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The possibility of different interpretations of the stochastic term (or calculi) in the overdamped 
Langevin equation for the motion of a particle in an inhomogeneous medium is often referred to as 
the "Ito-Stratonovich dilemma," although there is, in fact, a continuum of choices. We introduce two 
Monte Carlo (MC) simulation approaches for studying such systems, with both approaches giving 
the choice between different interpretations (in particular, Ito, Stratonovich, and "isothermal"). 
To demonstrate these approaches, we study the diffusion on a ID interval of a particle released 
at an interface (in the middle of the system) between two media where this particle has different 
diffusivities (for example, two fluids with different viscosities). We consider the properties of the 
particle distribution for reflecting boundary conditions at the ends of the ID interval. A discontinuity 
at the interface in the stationary-state particle distribution is found, except for the isothermal case, 
as expected. We also study the first-passage problem using absorbing boundary conditions. Good 
agreement is found when comparing the MC approaches against theoretical predictions as well 
as Brownian and Langevin dynamics simulations. Additionally, while this problem was chosen 
primarily to verify the algorithms, the results themselves turn out to be interesting — particularly 
when comparing across interpretations. For instance, we report that: 1) for some calculi, there 
can be more particles on the low-viscosity side at earlier times and then more particles on the high- 
viscosity side at later times; 2) there is no preference to end up on a particular wall for the Ito variant, 
but a bias towards the wall on the low-viscosity side in all other cases; 3) the mean first-passage time 
to the wall on the low-viscosity side grows as the viscosity on the high-viscosity side is increased, 
except for the isothermal case where it approaches a constant; 4) when the viscosity ratio is high, 
the first-passage-time distribution for the wall on the lower-viscosity side is much broader than for 
the other wall, with a power-law dependence of the former in a certain time interval whose exponent 
depends on the calculus; 5) the average portion of time the particle spends on a particular side can 
be very different from the probability to reach the wall on that side and depends significantly on 
how the averaging is done. 

PACS numbers: 02.70.Tt, 05.40. Jc, 05.10.Ln, 05.10.Gg 
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I. INTRODUCTION 

While deceptively simple to consider conceptually, the 
diffusion of a particle in an inhomogeneous medium is a 
subtle problem For a system consisting of a high 

viscosity and a low viscosity region, will a tracer particle 
be found more often on the "thicker" side, the "thinner" 
side, or with an equal probability at all locations? While 
an assumption of a Boltzmann distribution would suggest 
an equal probability at all locations, a simple physical ar- 
gument suggests that the particle may be more likely to 
be found on the high viscosity side. If the only inhomo- 
geneity in the system is the viscosity, then statistically, 
the dynamics are the same at all points with the excep- 
tion that time is rescaled on the high viscosity side where 
the particle moves slower. Hence the particle will spend 
a greater amount of time on that side and should be more 
likely to be found there. 

While this discussion can go back and forth, in truth 
each answer corresponds to a different formulation of the 
problem — both with their own mathematical and phys- 
ical assumptions. In essence, the question boils down 
to considering a particle at the viscosity interface: will 
there be an additional "push" at the interface to the 
low viscosity side that will compensate the "trapping" 



in the high viscosity region due to the slower dynamics 
there? Let us assume that the fluids can be considered 
implicitly, so their effect on the diffusing particle is re- 
duced to a random force and a friction force that the 
particle feels. We further assume that the interface is 
sharp and that there are no interfacial effects such as sur- 
face tension that would make the particle prefer or avoid 
the interface. Finally, we assume the particle motion is 
overdamped such that it can be described by the over- 
damped Langevin equation [f|. The situation where the 
assumption of equivalent dynamics (apart from the time 
scale) holds and particles accumulate on the high viscos- 
ity side corresponds to interpreting the stochastic term 
in the overdamped Langevin equation according to the 
rules of the Ito calculus Q. In other interpretations, of 
which the most well-known corresponds to Stratonovich 
calculus p} , but perhaps of most interest to physicists is 
the so-called isothermal formulation lUIllh the particle 
"feels" the interface and has a tendency to move to the 
low viscosity side. In particular, in the isothermal case 
(also referred to in the literature as Hanggi-Klimontovich 
or kinetic [l2j ) this bias at the interface cancels the trap- 
ping effect and, in the stationary state, the particle has 
an equal probability of being found at all locations. 

There have been several experimental explorations of 
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systems with nonuniform viscosity or, more generally, 
particle diffusivity [l|, Q ■ Most relevantly, van Milligen 
et al. [1] performed experiments in which the relative 
concentration of dye was monitored in a system com- 
prised of two viscosity regions created by adding gelatine 
to one half of it. In this setup, the dye was found to 
accumulate on the high viscosity side. As demonstrated 
in their paper, the concentration distributions they ob- 
served can be approximated by solutions of the Fokker- 
Planck equation, which corresponds to the Ito formula- 
tion (see Sec. HIl). On the other hand, in experiments by 
Lancon et al. [2], particles diffused between two nearly, 
but not perfectly, parallel walls, which, as the authors 
showed, is equivalent to diffusion with a variable diffu- 
sion coefficient. In contrast to Ref. while individ- 
ual particles were observed to drift to regions of higher 
diffusivity, there was no overall flux when starting with 
a uniform concentration of particles and hence the uni- 
form concentration was preserved, which is consistent 
with the isothermal case. On the theoretical side, the 
general conclusion has been that different formulations 
can be appropriate depending on the details of the sys- 
tem. For instance, in a simple model of diffusion in a 
slowly modulated periodic potential by Sokolov [12j , Ito, 
Stratonovich, and isothermal rules correspond to differ- 
ent variants of the model. 

Given the experimental and theoretical relevance of 
diffusion in inhomogeneous media, it is desirable to be 
able to model it computationally. One popular computa- 
tional approach to modeling the dynamics of particles in 
a fluid medium is molecular dynamics (MD) [13]. Two 
of the most common MD approaches assuming an im- 
plicit fluid are Brownian dynamics (BD) and Langevin 
dynamics (LD). The BD approach consists in solving nu- 
merically the overdamped Langevin equation in the Ito 
formulation. LD drops the overdamped assumption and 
adds an inertial term to the equation of motion that is 
then solved numerically. Interestingly, even in the limit 
of an infinitely large damping (but decreasing the time 
step so it is always much smaller than the corresponding 
relaxation time), LD is not reduced to BD and instead 
corresponds to the isothermal formulation of the over- 
damped Langevin equation. Thus the choice between 
the two methods is not arbitrary when an inhomogeneous 
system is considered. 

Monte Carlo (MC) methods are another frequently 
used group of computational approaches [l4j. Unlike 
MD, MC methods generally sacrifice at least some as- 
pects of the dynamics of the system to achieve a com- 
putational speedup compared to MD. How much is sac- 
rificed can vary. The popular Metropolis algorithm (l5j 
is designed solely to reproduce the equilibrium state and 
is generally unsuitable for studying the dynamics, even 
qualitatively. For instance, in the problem of inhomoge- 
neous diffusion considered here, a particle would diffuse 
equally rapidly regardless of the viscosity, unless care is 
taken to adjust the step length and/or the time incre- 
ment at every MC step (the time step) appropriately, 



which is not a part of the basic algorithm. However, 
there are also methods {dynamical MC) which do repro- 
duce the principal aspects of the dynamics and even the 
correct time scale of the process. This includes gener- 
ally faster, but sometimes less accurate and/or reliable 
methods where the particles are restricted to pre-defined 
points in space forming a lattice (lattice MC, or LMC), as 
well as off-lattice methods. Dynamical LMC methods in- 
clude kinetic Monte Carlo |16l4l8l | and methods designed 
to reproduce the correct dynamics in an arbitrarily strong 
uniform field [lMH. 

In this article, we present two dynamical LMC algo- 
rithms which can be used to generate results consistent 
with different formulations of the inhomogeneous diffu- 
sion problem. We consider a simple system consisting 
of a single random walker diffusing in a box which con- 
tains two separated fluids in equal proportions: a fluid 
of viscosity tjl on the left side and viscosity 77^ on the 
right side, similar to the experimental setup in Ref. 
Hence, there is a viscosity interface at the center of the 
box. Projecting the motion of the random walker onto 
the direction perpendicular to the interface reduces the 
problem to a ID one. The particular choice of the prob- 
lem with a sharp interface has been made for two (in a 
way, opposite) reasons. On the one hand, a priori it is the 
simplest setup in which to study diffusion in an inhomo- 
geneous medium; it lets us define simple quantities char- 
acterizing the process, such as the fraction of the time 
spent on a particular side of the interface. On the other 
hand, when it comes to actually solving the problem, it 
turns out there are difficulties specifically in the case of 
a sharp interface, as discussed in Sec. [TTJ so in that sense 
this problem is, in fact, more complicated than that of a 
gradual change in viscosity; our algorithms should gener- 
alize straightforwardly to the latter case. In addition, in 
future work we intend to apply the algorithms developed 
here to studying translocation of a polymer through a 
nanopore in a thin membrane, in which case having dif- 
ferent viscosities on the two sides of the membrane with 
a sharp boundary between them corresponding to the 
pore is natural. Overall, the system of a single random 
walker in a fluid with a single, sharp viscosity interface 
is chosen as a simple "toy problem" such that we can 
clearly demonstrate the validity of our LMC algorithms. 
However, as will be shown in the results section, many 
interesting features arise for the dynamics of a particle 
in such a system — particularly when contrasting results 
for different calculi. This manuscript thus not only in- 
troduces LMC approaches for simulating diffusion in a 
system containing regions of differing diffusivity, but it 
demonstrates the complexity which arises for even the 
simplest case and highlights the importance of choosing 
the calculus that is appropriate for the problem under 
study. 

In the next section, we briefly review the relevant parts 
of the theory of the overdamped Langevin equation for 
the motion of a Brownian particle and obtain the cor- 
responding equation for the particle concentration. We 
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then introduce Ito, Stratonovich, and isothermal calculi 
giving some examples of their physical relevance. We 
also point out some problems with using the overdamped 
Langevin equation in the case of a sharp interface. In the 
following two sections, we consider different simulation 
techniques that can be used for the problem of parti- 
cle diffusion. We first note that the well-known BD and 
LD methods correspond to Ito and isothermal calculi, re- 
spectively, and then introduce two MC algorithms. The 
details of the simulation protocol are outlined in Sec. V, 
after which in Sec. VI we present comparisons between 
different methods (BD, LD, MC) and results of the sim- 
ulations. We end the paper with a discussion of possible 
uses of the methods developed here and their future gen- 
eralizations. 



It is this limit that we will consider in what follows. How- 
ever, neglecting the inertial term can bring about subtle 
problems, as discussed below. 

Taking the case of free diffusion \U{x) = 0] we can 
rewrite Eq. §5§ as 



<2fcT 



dt 



R(t). 



(6) 



The particle diffusion coefficient (or diffusivity) D is de- 
fined by the expression 



(Ax 2 



2DAt, 



(7) 



where Ax is the particle displacement over time At. 
Based on Eq. ((HJ), it can be shown that 



II. THEORY 



A. Diffusion in a liquid 



We then obtain 



D 



kT 
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kT 
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(8) 



Diffusion of a pointlike particle along the x axis under 
the influence of a potential U(x) can be described by the 
Langevin equation 22] , 



dt 



(9) 



s dt dx w 



V2kTCR(t), (I) 



where m is the particle mass, k is the Boltzmann con- 
stant, T is the temperature, and R(t) is a random func- 
tion that satisfies (R(t)) = and (R(t)R(t r )) = 6{t - t'), 
where (. . .) denotes the ensemble average. The coefficient 
of the last term follows from the fluctuation-dissipation 
theorem. It is assumed that the drag force Fj (the first 
term on the right-hand side) is proportional to the par- 
ticle velocity v, 



F f = -C«, 



(2) 



where £ is the friction coefficient. For a particle moving 
in a fluid, £ is proportional to the fluid viscosity rj: 



(3) 



where K is constant for a given particle. For instance, for 
a spherical particle of radius 1Z, according to the Stokes' 
law (which assumes no-slip boundary conditions for the 
fluid), 



6ttK: 



(4) 



for other particle shapes and boundary conditions the 
expression for k will be different, but Eq. <j3j) is still valid. 

In many cases of interest, the left-hand side (or the 
inertial term) in Eq. ([1]) is small. If this term is ne- 
glected (the overdamped limit), we obtain the overdamped 
Langevin equation 



r d 



-±U(x) 
ax 



^2kT(R(t). 



(5) 



B. Inhomogeneous media 



The consideration below and in Section III CI is similar 
to Ref. [H, but somewhat more detailed. 

If we consider diffusion in a stationary inhomogeneous 
medium, D will be dependent on x such that 



—x = y/2D(x)R(t). 



(10) 



To proceed, we must determine how the multiplicative 
noise term in Eq. (|10l) is to be evaluated. For a deter- 
ministic differential equation 



dx 
It 



(11) 



where /(x, t) is a nonrandom and well-behaved function, 
we can write down the derivative on the left-hand side as 
the limit of a difference and then 



x(t + St) « x(t) + f(x(t), t)5t. 



(12) 



In the last term, instead of the value of the function / at 
time t, we could have taken its value at time t + St, or any 
linear combination of the two — this does not matter in 
the limit St — > 0, i.e., in this limit the equation 



x(t + St) 



x(t) + [(1 
+af(x(t-\ 



-a)f(x(t),t) 
St),t + St)]St 



(13) 



can be used with any a between zero and one. Varying 
a changes the right-hand side of Eq. (fT3|) by a negligi- 
ble amount 0(5t 2 ). However, for stochastic differential 
equations, like Eqs. © and (ITU1) . there are complications. 
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First, R(t) is not a well-behaved function. Instead of its 
instantaneous value, which is not well-defined, one needs 
to take its average over the interval St, which depends 
on St and is r/V&t, where r is a random number with 
(r) =0, (?' 2 ) = 1 and no correlation between different 
time intervals. Then, for Eq. (|9]), 



its average value 1, since the contribution of the fluctuat- 
ing part over the interval At is oc y/~N St — V AtSt, which 
is negligible when St 0. This gives 



x(t + St) &x(t) + rV2D8t. 



(14) 



Second, when D is variable [as in Eq. (jlQI) ]. there is an 
additional uncertainty, since the result depends on where 
D is evaluated, that is, on the value of a in 



x(t + St) « 



x(t) 



t)D(x(t)) + aD(x(t + St))]St, (15) 



and this dependence now remains even in the limit St — > 
0. Indeed, let us assume for now that D(x) is a smooth 
function. From now until the end of this section, we will 
use x(t + St) to denote the value calculated using Eq. 
(fl"5j) . rather than the actual value of the coordinate at 
time t + St, thus, the sign in Eq. (IT51) is replaced by 
the "=" sign. Expanding 

D(x(t + St)) - D(x(t)) = ^-[x(t + St) - x(t)} + 0(St) 

ax 

_ dD_ 

dx 
_ dD 

dx 
dD 



rV2[(l - a)D(x(t)) + aD(x(t + 5t))]St + 0{St) 
r^2D(x(t))St + 0(St 3 / 2 ) + 0(St) 



dx 



ry/2D{x{t))8t + 0(5t) 



(16) 



and using this in Eq. ([T5|) . we get 

x(t + St) - x(t) 



D(x(t)) + a(dD/dx)r v / 2D(x(t))St + O(St) 



St 



= ry/2D(x{t))St 



a(dD/dx)r y / 2D(x(t))St 
2D{x(t)) 



r^2D(x{t))5t + r 2 a(dD/dx)St + 0(St 3/2 ). 



0{St) 



x(t + St) « x(t) + ry/2D(x(t))St + u(dD/dx)St. (18) 

It is easy to show that the same result is obtained, if in 
Eq. (|T5|) the value of D is taken at a single point between 
the end points of the jump, as in Ref. [2|, that is, Eq. (fT5")) 
is replaced with 

x(t + St) = x(t) 



+r^2{D[x = (1 - a)x(t) + ax(t + 5t)]}St. (19) 

The last term in Eq. (|18[) is deterministic and produces 
drift of the particles in the direction of the diffusivity 
gradient known as "spurious flow" [23[ , "spurious drift" , 
or "noise- induced drift" [Ej. The drift velocity is 



V 



dD 

dx ' 



(20) 



it vanishes when a = 0. Note that the same drift can 
be produced by a suitably chosen external deterministic 
force, so it is always possible to change the value of a and 
still get the same dynamics by introducing an additional 
fictitious force 1 1211 . 



C. Particle fluxes and continuum diffusion 
equations 

Although the drift velocity of particles is zero for a = 
[Eq. (|2Ull ]. the net flux through a fixed point Xq does not 
vanish in this case, even without a particle concentration 
(or probability density) gradient. Indeed, assuming for 
simplicity that r = ±1 (with no influence on the final 
result), in a single jump the point with x — xq is crossed 
from the left by one half of all particles (namely, those 
with r = +1) with x between x\ and xq, where x\ obeys 
the equation 



(17) 



= xi + y / 2D(x 1 )St, 



(21) 



In Eqs. (fT6|) , (flTf and below, we have not specified ex- 
plicitly the point at which the derivative dD /dx is taken. 
It is easy to check that replacing the value of dD/dx at 
the point x(t) with that at the point x(t + St) or vice 
versa introduces an extra term of order St 3 / 2 in the last 
line of Eq. (fT7| that can be neglected, so the exact point 
at which the derivative is taken does not matter. On the 
other hand, both the first and the second term in the last 
line of Eq. (|17p are important in the limit St — > 0, even 
though they are of different orders in St, because (r) = 0, 
but (r 2 ) = 1^0. Indeed, over a fixed time interval At 
such that At = NSt, the contr ibuti on of the first term as 
a function of N and St is cx V NSt — VAt, while that of 
the second term is cx NSt — At, and both terms survive 
when St — > 0. In the second term, we can replace r 2 with 



xq — x\ ss \/2[D{xq) — (dD/dx)(xo — xi)]St 

w ^2D{x a )5t (l - —L-^(x - Xl )) , (22) 



V 2D(x ) dx 



with the solution 



x - Xl 



dD/dx 
y/2D(x ) 



Vsi 



« v /2D(x )St- (dD/dx)St. (23) 
The number of such particles is 

~ (y/2D(x )8t - (dD/dx)5t) p, (24) 



n+ 



■5 



where p is the (uniform) particle concentration. Likewise, 
the same point xq is crossed from the right by one half of 
all particles (those with r = — 1) with x between Xq and 
X2, with 



x 2 - x « y/2D{x a )St + (dD/dx)St, (25) 
and the number of such particles is 

n_ fa ~ (' A /2 J D(a;o)^ + (dD/dx)5t) p. (26) 
The net flux is thus 



(n+ - n_)/£f = -p^. 



(27) 



When a ^ 0, the flux due to the drift, pV, with V given 
by Eq. ([20]), has to be added to Eq. ([27]) . which gives the 
total flux 



. dD 

J=-(l-a)p—. 

ax 



(28) 



This flux is nonzero when the drift velocity is zero at 
a = 0, but it vanishes when a = 1. When a concentration 
gradient is present [p = p(x,t)], this adds the usual flux 
given by Fick's first law, —D 



dp. 

dx ' 



t n \ dD dp 
J = - {l - a)p ^- D d-x 



(29) 



In two special cases this expression is simplified: for a — 
0, 



J = — 



d(Dp) 
dx 



while for a = 1, 



ax 



(30) 



(31) 



which conforms to normal Fickian diffusion rules. The 
rate of the concentration change is 

dp dJ , . d ( dD\ d ( »p\ , nn . 

i = -^x = ^- a) Tx{ p ^) + d-x\ D Tx)- (32) 

In particular, for a = 0, 

dp _ d 2 (Dp) 



Of 



dx 2 



(33) 



which is the standard form of the Fokker-Planck equation 
without the drift term 23]; for a = 1, 



dp = d_ /dp 
dt dx \ dx 



(34) 



which is Fick's second law. 

The stationary state can be obtained by setting 
dp/dt = in Eq. ([32]). This condition for no net change 
in the particle concentration is then 

- J = (1 - a)p^ + D°f- = const. (35) 
dx dx 



If there is no net flux, as, for example, when the system is 
confined between two reflecting walls, then the constant 
in Eq. (1351) is zero and 



(1 - a)pdD = -Ddp, 



which gives 



pD 1 ' 



const. 



(36) 



(37) 



For a = 1, the particles can be found with equal proba- 
bility anywhere (p = const), but for a < 1 they are more 
likely to be found in regions of lower diffusivity (or higher 
viscosity) . 



D. Ito, Stratonovich, and isothermal calculi 

Even though the actual solution of Eq. ([TU]) depends 
on the value of a in its discretized variant [Eq. ([15]) ]. its 
formal solution can always be written as 



x(t) = x(0) + / y / 2D{x(t>))R{t')dt'. 



(38) 



Different a thus correspond to different interpretations of 
the formal integral in Eq. (I3"51) , or different calculi. Three 
cases have received special attention. 

The case a — corresponds to Ito calculus @ . In this 
case, the stationary distribution is given by 



p{x)D{x) — const. 



(39) 



The primary attractive quality of this formulation is 
that its implementation is "non-predictive": According 
to Eq. (fT"5]) . the value of the diffusivity D is taken at 
the initial point of the jump only, so, as opposed to all 
other cases, one does not need to look ahead in time 
to evaluate the spatially dependent term. This property 
makes the Ito formulation a frequent choice for mathe- 
matical problems — including financial modeling. But it 
is also physically relevant: For instance, the overdamped 
limit of the Langevin equation (JTJ corresponds to the Ito 
calculus in the idealized situation where the temperature 
varies in space, but the friction coefficient is temperature- 
independent and thus the same everywhere 24j. This is 
an inherently non-equilibrium situation from the thermo- 
dynamic point of view and so the fact that the station- 
ary distribution does not coincide with the equilibrium 
Boltzmann distribution is not surprising. Another ex- 
ample would be a particle in a modulated quasiperiodic 
potential where the tops of the barriers are all at the same 
level, but the well depths are position-dependent [12| . In 
this case, the change in diffusivity is accompanied by the 
change in the free energy and the Boltzmann distribution 
is obeyed. 

The case a = 1/2 corresponds to the Stratonovich cal- 
culus Q- In this case, the value of D at the midpoint of 
the jump (or the average of the values at the endpoints) 
is taken. Stratonovich calculus conforms to the "normal" 
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rules of calculus. For example, if x obeys Eq. (fTUj) . then 
one would expect 

^-x 2 = 2x^- = 2xJ2D(x)R(t), (40) 
at at 

where the right-hand side is interpreted according to 
the rules of the given calculus. This is indeed true for 
Stratonovich, but not for other calculi. The Stratonovich 
calculus arises when considering the Langevin equation 
(HJ) with the noise R(t) having a finite and space- and 
time-independent correlation time, in the limit when 
both this correlation time and the relaxation time m/C 
go to zero, but the former is kept much larger than the 
latter 25]. In fact, even for the overdamped Langevin 
equation (fTUj) (corresponding to the relaxation time equal 
to zero) with a finite noise correlation time the issue of in- 
terpretation of the noise term does not arise and the limit 
when the correlation time goes to zero corresponds unam- 
biguously to the Stratonovich calculus 0, [23| . This situ- 
ation violates the fluctuation-dissipation theorem which 
requires that a finite noise correlation time be accompa- 
nied by a memory in the friction force, thus giving rise to 
the generalized Langevin equation [22J instead of the or- 
dinary one. So the fact that the density in the stationary 
state is not uniform, but rather obeys 

pVD = const, (41) 

is natural. Sokolov [l2[ has also proposed a variant of 
his modulated quasiperiodic potential model that cor- 
responds to the Stratonovich calculus; again, the free 
energy in that model is position-dependent. Since Ito 
and Stratonovich are the two calculi most often dis- 
cussed in the literature, the choice between different cal- 
culi has often been referred to as the "Ito-Stratonovich 
dilemma" [H. 

But besides these two cases, the case a = 1 is also of 
special interest, because in that case Eq. ([3"T]) becomes 

p = const, (42) 

which corresponds to the equilibrium Boltzmann distri- 
bution in the "purest" case when it is only the diffusiv- 
ity that varies in space, but the particle free energy and 
the temperature are the same everywhere. There is no 
universally accepted name for the a = 1 calculus, with 
names such as Hanggi-Klimontovich , 12j , kinetic [ill, E3 
and isothermal @ appearing in the literature. We will 
use the latter in what follows, since it emphasizes the 
fact that this calculus corresponds to the case when the 
system is kept at a constant temperature, as opposed to 
the Ito case which can be achieved by varying the tem- 
perature across the system (still assuming the situation 
where the free energy would be constant for T = const). 

E. Discontinuous D(x) 

So far, we have assumed that D(x) varies smoothly in 
space. First of all, this allowed us to do the expansion in 



Eq. (|16[) . Even more importantly, we note that if there 
is a jump in the diffusivity, then even the solution of 
Eq. ([TO)) is not well-defined for a ^ 0. Consider a — 1 
and let us assume that 

D w = { %. I > I: <«> 

with Dl > Dfj. Let the particle coordinate at time 
t be xo < and suppose that the time step 5t and 
the random number r chosen at the following step are 
such that x + y/2D L rSt > 0, but x + y/2D R rSt < 0. 
Note that the jump size depends on the final position 
at time t + St. If we assume that the particle ends up 
to the right of the interface, then the jump size should 
be y/2DnrSt, but then the particle coordinate after the 
jump is xq + \/2Dnrdt < 0, so we come to a contradic- 
tion. Conversely, if we assume that the particle ends up 
on the left side, then the jump size should be \/2DirSt, 
but Xq + \J2Di,rSt > 0, so there is again a contradiction. 
Thus at this point the particle behavior is undefined. On 
the other hand, for xq > it is possible that both as- 
sumptions lead to valid results, so there is ambiguity. 
Given that in principle |xo| can be arbitrarily close to 
zero, this will always be a potential problem, no matter 
how small St is. In fact, this problem exists whenever the 
jump size depends not just on the initial position, that is, 
for any a ^ 0. Thus, in order to be able to use Eq. (fTUj) 
[discretizing it as in Eq. (|15p ] for calculi other than Ito 
(and even for the very concept of different calculi to be 
meaningful), we need to assume that the function D{x) 
is smooth everywhere. In that case, what we mean by 
a "sharp" interface and a "jump" in diffusivity is that 
actually the diffusivity is smooth and thus the interface 
between two liquids has a finite width, but the interface is 
narrow compared to all other length scales (in our simple 
problem, there is only one such length scale, the size of 
the system). Physically, the smoothness of D means that 
it changes little on the length scale of the typical inter- 
atomic distance. This is very likely to be true in reality, 
and even if the liquids are completely immiscible and the 
interface is atomically sharp, but the diffusing particle is 
much larger than the interatomic distance, the proper- 
ties of the system are effectively averaged out over the 
particle size and the condition is satisfied. In this case, 
Eq. (ITU1) is still meaningful for all calculi and in princi- 
ple can be solved numerically by choosing a very small 
time step when iterating Eq. (|15[) . so that typical parti- 
cle jumps in the simulation are smaller than the width of 
the interface. However, this would be a waste of compu- 
tational resources; but choosing a larger step would lead 
to the indeterminacy problem described above. This is 
where LMC methods are helpful, as we will see below. 

In the remainder of this paper, we will concentrate on 
this situation of a sharp, but not infinitely sharp inter- 
face. For simplicity, we will still use notation that may 
seem to imply that the interface is infinitely sharp, but 
will keep in mind that in reality it is not. So, for instance, 
p(+0, t) denotes the particle concentration immediately 
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to the right of the interface, but still far enough from 
it that the diffusivity has already reached its asymptotic 
value Dr; likewise, D(+0) = Dr. Bearing this in mind, 
Eq. ([57)1 for the stationary distribution and its particular 
cases, Eqs. (j3"9")l . (HI]), and (H2"1) . should remain valid in 
our situation. Moreover, we can show that these rela- 
tions are valid at all times, even in the transient state, in 
the vicinity of the interface, that is, at all times 



p(-0,t)[£>(-0)] 1 - Q = P (+0,t)[D(+0)Y 



(44) 



Indeed, with the exception of the case a = 1, the first 
term on the right-hand side of Eq. (|2T))) for the particle 
flux is large at the interface and in fact infinite in the limit 
of an infinitely sharp interface. However, the flux itself, 
even though it can be nonzero when the stationary state 
is not reached, is expected to be finite. Therefore, the 
second term on the right-hand side of Eq. (l29l) should be 
large (infinite in the limit) as well, and the flux J can be 
neglected compared to either of these terms. Dropping J 
in Eq. (j2"9l produces the stationary-state equation (f5Fj|). 
from which Eq. (j44j) follows immediately. In the spe- 
cial case a = 1, the first term on the right-hand side of 
Eq. (|2T)1) is zero and therefore the second term is finite, 
thus there is no jump in p, which is likewise consistent 
with Eq. (|44p . Condition (|44|) describes the concentra- 
tion jump across the interface that forms immediately 
and always exists, if a ^ 1. The full solution of the dif- 
fusion problem for the system consisting of two uniform 
regions separated by a "sharp" interface (where "sharp" 
is understood as above) can be obtained by solving two 
uniform diffusion problems in the two regions, 



dp(x,t) 
dt 



D 



L,R- 



d 2 P (x,t) 
dx 2 



(45) 



with two additional matching conditions at the interface: 
Eq. (|44[) and the flux continuity condition, 



D, 



dp{x,t) 



dx 



= D f 



dp{x,t) 



x=-0 



dx 



(46) 



a:=+0 



In addition, of course, there have to be boundary condi- 
tions at the walls, for example, the no-flux condition 



dp(x, t) 



dx 

for a reflecting boundary or 







p{x b ,t) = 



(47) 



(48) 



for an absorbing boundary, where Xb is the coordinate of 
the boundary. 



III. MOLECULAR DYNAMICS SIMULATION 
APPROACHES 

In this section and the next one, the various simula- 
tion approaches used to explore diffusion in an inhomo- 
geneous medium are presented. We start by describing 



two well-known molecular dynamics (MD) approaches, 
namely, Langevin dynamics (LD) and Brownian dynam- 
ics (BD). 

MD is a simulation approach based on solving equa- 
tions of motions for particles comprising the system. In 
the most straightforward, but very computationally in- 
tensive approach, both the probe particle(s) whose dif- 
fusion is studied and the fluid particles are included ex- 
plicitly and Newton's Second Law equations 



nriiXi — Fi 



(49) 



are solved numerically, where rrii and the mass and 

the position of the ith particle and Fi is the sum of all 
forces acting on it. The computational expense can be 
reduced by considering the fluid particles implicitly, as is 
done in the LD and BD approaches. 



A. Langevin Dynamics 

In the LD approach ((26[, Sec. 14.4), only the equation 
of motion of the diffusing particle is considered explicitly 
and the action of the fluid particles on it is replaced by 
the friction force and a random force, and the resulting 
equation of motion is given in ID by Eq. ([1]) with U = 0: 



mx = 



(50) 



Dividing by £ and assuming that £ is a function of x, we 
obtain 



or 



m 



mD(x) 




= -x + y/2D(x)R(t). 



(51) 



(52) 



Unlike in the case of the first-order Eq. ((TO]) , there is no 
ambiguity due to different interpretations of the random 
term even when £ (or D) is space-dependent, with all in- 
terpretations and all "correct" numerical methods giving 
the same solution in the limit when the time step St — » 0. 
This is because, unlike in the case of Eq. (fT0|l . the parti- 
cle velocity does not diverge, thus when the equation is 
discretized, the displacement during a single step is much 
smaller [0(5t) instead of 0(8t 1 / 2 )], so in the limit St — > 
it no longer matters whether the position at the begin- 
ningor at the end of the step is used when calculating 
D |24l|. In practice, some variation of the Verlet method 
([261, Sec. 13.4) is usually used to solve Eq. (|52|) numer- 
ically. By construction, Eq. (|52p obeys the fluctuation- 
dissipation theorem, thus in the case T = const the sta- 
tionary distribution should coincide with the Boltzmann 
distribution p — const. Since this fact does not depend 
on the value of the parameters entering the equation, it 
should remain valid in the limit m — > when the incrtial 
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term vanishes, and thus this overdamped limit should 
correspond to the isothermal interpretation of Eq. (fT0|) . 
In fact, this was proved in Ref. [27| . However, in order 
to reproduce this result numerically, it is important to 
retain the inertial term decreasing its coefficient, rather 
than simply dropping it, and to make sure that the time 
step is always much smaller than the relaxation time, 
m/(. Remarkably, unlike Eq. (dUJ) with a/0, Eq. (j52]l 
can still be solved straightforwardly in the case of a sharp 
interface, as we will see in Sec. VI. 



B. Brownian Dynamics 

The motion of microscopic particles in a liquid is nearly 
always overdamped, with the time of relaxation to the 
terminal velocity, m/£, usually much smaller than any 
other relevant time scale. As mentioned above, the time 
step of the LD algorithm should be much smaller than 
even this extremely small time. Even if the particle mass 
m is increased artificially to increase the relaxation time 
and thus the maximum allowed time step, this time step 
may still end up being too small and as a result, the sim- 
ulation will be inefficient. A seemingly straightforward 
approach is eliminating the left-hand side of Eq. (|52"j) (the 
inertial term) altogether and solving the resulting first- 
order equation. This equation coincides with Eq. (ITUl) . 
so the discussion of the ambiguity of the stochastic term 
in Sec. [TT| applies. Normally in MD simulations, to solve 
Eq. (ITUl) numerically the simplest Euler algorithm (of- 
ten referred to as Euler-Maruyama in the stochastic case 
fig! ) is used, in which, in order to obtain the particle po- 
sition at time t + St, the right-hand side of the equation 
is evaluated at time t. This approach, called Brownian 
dynamics, obviously corresponds to the Ito interpreta- 
tion of the stochastic differential equation (|10p . Thus, 
surprisingly, LD (even in the overdamped limit) and BD 
will give different results when the friction coefficient C, 
varies in space, as they correspond to different calculi. 
Just as LD, BD can be used in the case of a sharp inter- 
face without any complications, since these complications 
only arise for a/0 when solving Eq. (fT0|) . Some diffi- 
culties arise when the random number r is chosen to be 
discrete (e.g., ±1), as described in Sec. II V Cl so using r 
with a continuous distribution is preferable. 



IV. MONTE CARLO APPROACHES 

The main goal of this paper is to develop and test LMC 
algorithms for simulating a random walker at a viscosity 
interface corresponding to different calculi [a in Eq. (p~5|) ] . 
including the popular Ito, Stratonovich, and isothermal 
calculi. There are several reasons to develop such algo- 
rithms. First, LMC algorithms are popular for studying 
diffusion problems |29l445j and in some of these problems 
studying the case when the effective viscosity is inhomo- 
gencous in space (e.g., different degrees of crowding in 



different parts of the system) can be of interest (see, e.g., 
Ref. (30[). Second, LMC algorithms can offer advantages 
compared to MD (or, in general, solving numerically the 
Langevin equation) specifically in the case of a sharp in- 
terface. As we have seen, sharp interfaces are algorithmi- 
cally problematic when solving the overdamped Langevin 
equation, except in the case a = corresponding to BD. 
While the a = 1 (isothermal) case can be treated by LD, 
this may require a very small time step for the solution 
to be reliable. All other cases are in principle equiva- 
lent to Ito plus an additional drift force, but that drift 
force would be infinite at the interface and zero elsewhere, 
which would cause obvious numerical problems. LMC al- 
gorithms, on the other hand, are straightforward in all of 
these cases, as we will demonstrate. 

Two different approaches will be taken. In the first 
method, changes in viscosity values are introduced by al- 
tering the probability p s that a particle does not jump 
during a MC jump attempt: the particle has a higher 
chance of "staying put" in the high viscosity region. In 
this approach, the jump lengths on the low and high 
viscosity sides are the same. Conversely, in the second 
method, the change in viscosity across the interface is in- 
cluded by altering the jump length (i.e., the mesh step of 
the lattice): The particle jumps shorter distances when 
the viscosity is higher. Here, the probability of not jump- 
ing is constant (except perhaps for sites adjacent to the 
interface). In both methods the time step is the same 
in both regions and stays constant throughout the sim- 
ulation. This is essential for so-called numerically exact 
variants of LMC methods [HI, M, M, IH B H S3 and 
is also very helpful when several particles are simulated 
at once. 

Both LMC methods are derived based on three require- 
ments: (I) the diffusion should be unbiased away from 
the interface, (2) the diffusion rates should be correct in 
both regions, and (3) the probabilities of the moves at 
the interface should ensure that the stationary distribu- 
tion given by Eq. (|37|) is reproduced. Only the last of 
these requirements depends on the calculus (the value of 
a), and so this value will only affect the rules of the al- 
gorithms at the interface, that is, in the only place where 
the viscosity gradient is present, as expected. 

When deriving both algorithms, the following consid- 
eration is used. Suppose the particle moves on a lattice 
with a constant mesh step a and the time step is t. At 
every successful step (i.e., ignoring the time steps where 
the particle stays put), the particle moves left or right 
with an equal probability. The mean-square displace- 
ment after N successful steps is 

(Ax 2 ) = Na 2 . (53) 

Since in time At there are At/r total steps and thus on 
average At(l — p s )/r successful steps, 

{Ax2) = (1^^ (M) 
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and comparing this to Eq. 0, 

(1-PsW 



D = 



2t 



(55) 



A. MC I: Variable p s , fixed jump length a 

1. The basic algorithm 

In the first Monte Carlo approach, all "successful" 
jumps have the same length and the dependence on vis- 
cosity is introduced by altering the probability of remain- 
ing still during a time step (Fig. [TJ . 



Both 



pi 



Isothermal 



Ito 



ft 



A 



FIG. 1: A schematic of the MC I algorithm. Different shades 
of gray depict regions of different viscosity, with the darker 
region being more viscous. Vertical lines correspond to the 
locations of ID lattice sites. The distance between sites is the 
same in both regions. Arrows show transitions between sites: 
the thicker the arrow, the higher the transition probability. 
The notation for the probabilities corresponds to that used in 
the text. The top part illustrates the behavior of a particle 
away from the interface, which is independent of the calculus, 
while the middle and bottom parts are for a particle at the 
interface site in the two particular cases of the isothermal and 
Ito calculi, respectively. 

Given the diffusion coefficient D and the time step r, 
from Eq. (|55|) the probabilities to move right and left 
should be 



P+ = P- 



1-Ps 



2 a 
and the probability to stay put, 



Dt 

2 ' 



Ps 



1 



2Dt 



(56) 



(57) 



Obviously, r and a need to be chosen so that all of these 
probabilities are between zero and one for the diffusivitics 
both to the left (D = D L ) and to the right (D = D R ) 
of the interface. Since D is inversely proportional to the 
fluid viscosity rj [see Eq. ((5J], it follows from Eq. (|56|) 
that p + and p- should be inversely proportional to rj, or, 
introducing the viscosity n\ at which these probabilities 



are equal to 1/2 and correspondingly p s = 0, 



L,R L,R 

Pi+ = Pi- = 



vl 



L.R 



Pis 



1 



_Vo_ 

Vl,r 



(58) 
(59) 



where sub- and superscripts L and R denote the left and 
right side of the viscosity interface, respectively, and "I" 
shows that these parameters correspond to the MC I al- 
gorithm. The same probabilities were used previously 
in Refs. [48|, |49(. The parameter n\ should not be larger 
than the smaller of the two viscosities t]l and r] R to ensure 
that both and pf s are non-negative. In fact, even the 
situation where r]\ is equal to or only slightly smaller than 
either t)l or rjn should be avoided, since some simulation 
artifacts, such as spurious oscillations both in time and 
in space, are possible when the waiting time is zero |50| . 
or, more generally, the variance of the time between suc- 
cessful steps is small. On the other hand, from Eqs. © 
and (j56"j) . 



j _ 2kTn 



(60) 



so this limitation on the allowed values of rj\ restricts pos- 
sible choices of tj and oi. The fact that the probabilities 
of the moves are different in the two parts of the system 
is depicted schematically in the top part of Fig. [1] 

Equations {55| and (|5"9")) should be used for lattice sites 
away from the interface. They ensure that the first two of 
the three requirements mentioned above (no bias and the 
correct diffusivities) are satisfied. The third condition, 
the correct ratio of the particle concentrations on the two 
sides of the interface, should be achieved by choosing a 
special set of probabilities for sites near the interface. It 
is convenient to choose the lattice so that the interface is 
on a lattice site. Let us denote this site 0, and the sites 
immediately to the left and to the right of it, —1 and +1, 
respectively. Since we have a single condition to satisfy 
for several probabilities, we can make an arbitrary choice 
that the probabilities of moves out of sites —1 and +1 
are still given by Eqs. ([58]) and (|59|) and that only the 
two probabilities for the moves out of site (plus the rest 
probability at site 0) are special. Based on Eq. (|57|) . the 
condition that needs to be satisfied is that 



n_i 



Vr 



(61) 



where rii is the average number of particles at site i in 
the stationary state. Since from detailed balance, 



and 



ni 
n 



n 
n_i 



Pl+ 

Pi- 



Pt\ 

pi 



(62) 



(63) 
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where pj + and pj_ are the probabilities of the moves from 
site to the right and to the left, respectively, we get 



n-i 



or 



Pi- Pi- 



P°U 



pj+m 

Pi- VL 



VR 
VL 



Vr 

VL 



(64) 



(65) 



This means that for tjr ^ t]l and a ^ there is a bias 
at the interface {p\_ ^ Pi + )- 

Equation (|65|) defines the ratio but not each 

of these probabilities separately. In the isothermal case 
(a = 1), p = const in the stationary state and thus n\ = 
ri_i; it is reasonable to require that uq is the same as 
n_i and ni, and this gives unambiguously 



rf _( = 1 , = pf+ = A, 



The probability of staying put is then 



Pi 



where the reduced viscosity 



^VlVr 
Vl + 



(66) 
(67) 

(68) 
(69) 



is the harmonic mean of the two viscosities. Since it 
makes sense for the rate of moving away from the inter- 
face to always be related to some sort of average between 
the two viscosities, we can require, somewhat arbitrarily, 
that the probability of staying put be given by Eq. (|6"8")l 
for all values of a and use this as an additional condition 
together with Eq. (|6"51 to obtain 



o = Vo Vr 

v*vl + v a R 

n _ Vo Vl 
Pi+ 



V VI + Vr 



(70) 
(71) 



In fact, we have checked for the Ito calculus that the 
dynamics are insensitive to the choice of p® s provided 
that pj + = , which in this case follows directly from 
Eq. (|65l) . Results essentially indistinguishable from those 
generated by the above formulation were obtained in sim- 
ulations with j»j s — and Pi ± = 0.5. 

Equations (|70|) and (|7Tj) have simple interpretations in 
two cases. In the isothermal case (a = 1), pj_ = an d 
Pj + = pf±. In other words, the probability of the move 
to the left from the interface site (the move that is within 
the left region) is the same as the probability of any other 
move (from a non-interface site) within the left region, 



and the same is true for the right region. This property 
is illustrated in the middle part of Fig. [TJ It ensures 
explicitly the absense of a net flux through any surface 
when the particle concentration is the same everywhere, 
the key property of the a = 1 case mentioned in Sec. [Til 
since for every surface drawn between two sites there will 
be as many particles moving from left to right as from 
right to left. In the Ito case (a — 0), the probabilities of 
moves to the left and to the right from the interface site 
are equal (the bottom part of Fig. [l}. This is consistent 
with the description of the Ito formulation in Sees. U and 
HJas the case where any bias at the interface is absent. 



2. Variants with all move attempts successful 

When the ratio of the two viscosities, tjl and 77^, is 
very large or very small, then, since t]q should be chosen 
at least as small as the smaller of the two, the proba- 
bility of moving in the more viscous region will be very 
low, according to Eq. (1551) . This makes the algorithm 
very inefficient in this case, since many random numbers 
need to be generated per single successful move. How- 
ever, one may note that the number of steps until the 
next successful move is exponentially distributed, with 
the probability that this number is n being 



Pin = (1 -pis)Pi~ 



(72) 



where p\ s is either p^ s , p^ s , or pj s , depending on the lo- 
cation of the particle. Instead of waiting for many steps 
until a successful move, one can move at every step, but 
advance the clock by m\ instead of tj, generating n from 
the distribution (|72[) . The resulting procedure is exactly 
equivalent to the original one. Strictly speaking, in this 
procedure the time step is variable, which in some cases 
is undesirable, as mentioned above; however, it is still al- 
ways a multiple of the basic step tj . It is also possible to 
always advance the clock by the same amount ti w , equal 
to the average waiting time between successful steps, 



n 77 

1 - Pis Vo 



(73) 



where r\ is either 77^, 77^, or 77*. While not strictly equiva- 
lent to the algorithm given by Eqs. (|5"8"]l and (|5"9"]l . such an 
algorithm should ultimately give the same results in the 
limit a, r — > 0. The advantage is that the need to draw 
random numbers from the distribution (|72[) is avoided. 
However, the time step will no longer be a multiple of 
t. A potentially more serious problem is simulation arti- 
facts that may appear when the waiting time is constant, 
as mentioned above. Regardless of the utility of this al- 
gorithm for actual simulations, it is useful in the analysis 
of the diffusion problem, as we discuss below. 

In both variants of the algorithm (with random and de- 
terministic clock advancement), at each step the moves 
to the left and to the right are equiprobable, except when 
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the particle is right at the viscosity interface the proba- 
bilities are 



P?- 



Vr 



(74) 
(75) 



These probabilities are such that their ratio p'i + /p'i_ is 
the same as Pj + /pj_, but their sum is unity, so the rest 
probability is zero. 

Note that in the algorithm presented here the particle 
moves to a neighboring site at every step, that is, the 
particle motion is simply an ordinary random walk on 
a ID lattice, and it is unbiased everywhere, except per- 
haps at the interface. (We ignore for the moment the 
fact that each step takes a different amount of time in 
different regions.) In fact, in the special case of the Ito 
calculus the random walk is unbiased at the interface as 
well, since in that case Eqs. (|T4"|) and ([751 give p'®_ = p'® + . 
In general, the trajectory of the particle can be thought 
of as consisting of segments of an unbiased random walk 
terminating at the interface; every time the interface is 
reached, the particle chooses to continue on the left side 
with the probability p'®_ or on the right side with the 
probability p'i + . Each such choice is independent of all 
of the previous choices and of whether the particle has 
reached the interface from the left or from the right. We 
will use this consideration in Sec. VIB and in Appendix 
B to derive some properties of the system. 



B. MC II: Fixed p s , variable jump length 

In the second MC approach, the probability of remain- 
ing still, pus, is independent of the viscosity and instead 
the length of a successful jump is dependent on the vis- 
cosity. To keep the particle on-lattice, we introduce a 
mesh with a different step in the two regions (Fig. [2]). 
From Eqs. (|55|) and ([5]), when pu s and Tn are fixed, the 
mesh step a cx \f~D cx 1/ \/rj, so the mesh steps on the left 
and right sides are, respectively, 



ano 




where ano and ^J 1 are constants and 

2kTm 



n 11 



KCliloil-PUs) 



(76) 
(77) 

(78) 



While it is possible to put the interface at one of the 
lattice sites, as in MC I, we now place it between two sites 
to illustrate the possibility of this choice. The distance 
between the two sites closest to the interface (one on each 



side) is chosen to be the average of the mesh steps on the 
two sides, ajj = (afj + a/j)/2, although the exact value 
does not matter in the limit afj,a^ — > 0. 
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FIG. 2: A schematic of the MC II algorithm. Different shades 
of gray depict regions of different viscosity, with the darker 
region being more viscous. Vertical lines correspond to the 
locations of ID lattice sites. The distance between sites is 
smaller in the region of higher viscosity. Arrows show tran- 
sitions between sites: the thicker the arrow, the higher the 
transition probability. Notation for the probabilities corre- 
sponds to that used in the text. The top part illustrates the 
behavior of a particle away from the interface, which is in- 
dependent of the calculus, while the rest of the figure is for 
particles on sites adjacent to the interface on both sides in 
the cases of the isothermal, Ito, and Stratonovich calculi. 

For all sites not adjacent to the interface, on both sides, 
the probabilities of the moves are the same, with 



L,R. L,R 
Pu+ =Pu- 



(l-piis)/2 



(79) 



and the probability of remaining still, pu s , is the same at 
all these sites, as well. The choice of pu s is in principle 
arbitrary. While pu s — seems the most efficient choice 
at first glance, since the particle would move at every 
step of the algorithm, this, as mentioned, can introduce 
artifacts and, in fact, we have recently shown that 
piis = 2/3 is optimal from the point of view of accuracy 
of the algorithm and may even be more efficient than 
Piis = 0, since a coarser mesh can be used to achieve the 
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same accuracy. For this reason, we choose p\\ s = 2/3 in 
this work. 

As in MC I, the probabilities of the moves need to be 
modified for sites adjacent to the interface. A reasonable 
choice is to allow the probabilities of the two moves cross- 
ing the interface to change but keep the probabilities of 
the two other moves leaving these sites that do not cross 
the interface the same as given by Eq. (|75|) . We use the 
following notation for the probabilities of the moves that 
are allowed to change: denotes the probability of 

the move to the right from the site left of the boundary 
and is the probability of the move to the left from 
the site right of the boundary. This notation is similar 
to that we used for MC I, but it has a somewhat differ- 
ent meaning: in MC I these were the probabilities of the 
moves from the same interface site, but here these moves 
cross the interface starting from different sites. 

The detailed balance condition is 



n L 



pJ^ 

P°u- 



(80) 



where ul (jir) is the average number of particles per site 
to the left (right) of the interface in the stationary state. 
To relate this condition to the stationary distribution of 
Eq. (l37l) . we need to keep in mind that the density of 
sites (their number per unit length) is different in the 
two regions. The particle concentration p is proportional 
to both the number of particles per site and the site den- 
sity; the latter is inversely proportional to the mesh step. 
Therefore, using Eqs. (175)1 and ([77]) for the mesh steps, 



PR 

PL 




or 



pU 

P°n- 



VL 



1/2-a 



(81) 



(82) 



Again, as in MC I [cf. Eq. (|65p ]. we have one equation for 
two probabilities. As an additional condition, we choose 



max(p n+ ,p n _) =Pri± , 
which, together with Eq. (|52"j) . gives 



o l.r ■ 

Pn+ = Pii± min 



PlI- 



L,R- 
Pu± mm 



1/2-a 



1/2-a 



(83) 

(84) 
(85) 



This is a simple choice that guarantees that both pf I+ 
and are always between zero and one for any value 
of a and any choice of pu s . Note that always either 
Pi 1+ = Pi{± or = Pif±, so only one of these probabil- 
ities (at most) gets modified compared to the bulk value. 
The corresponding probability of staying put is changed 



accordingly. In the special case a = 1/2 (Stratonovich 
calculus), p1 1+ = = Pii'±, so all probabilities of all 
moves are the same (Fig. [2]). While this means that 
for a — 1/2 there is no flux across the interface when 
nR = til, this is not the same as pr — pl (a uniform 
particle distribution) . 

We note that a similar, but off-lattice approach was 
successfully used recently by one of the authors [5l| to 
treat a problem corresponding to the a = 1 case. 



C. Comparison between MC II and BD 

It is instructive to compare the MC II algorithm de- 
scribed here in the Ito case (a = 0) to the BD algorithm, 
a single step of which is given by Eq. (|15p with a = 0, 
specifically its variant with the random number r = ±1. 
Similarities are obvious: in both cases the time step is 
fixed and the jump size depends on what region the par- 
ticle is in, being proportional to \ D (or inversely propor- 
tional to y/rj). However, what happens at the interface 
is different for these two algorithms. In MC II, in order 
to stay on-lattice, the lengths of the jumps across the 
interface are chosen to be the same in both directions, 
but in BD the jump length depends on what region the 
jump originates in and thus is different for the left-to- 
right jumps and the right-to-left jumps. This difference 
is crucial and makes it possible to achieve the correct par- 
ticle concentration jump across the interface in the sta- 
tionary state while keeping the probabilities of all jumps 
the same (effectively equal to 1 /2, since there is no bias in 
BD and no staying at rest), whereas in MC II obtaining 
this concentration jump requires setting p^_ ^ p_. This 
is actually a rather subtle issue. Indeed, in the stationary 
state, the left-to-right and the right-to-left fluxes across 
the interface should be equal. But then one could argue 
for the BD algorithm that since the flux is equal to the 
product of the particle concentration and the jump size 
and the latter is proportional to VT), then the stationary 
concentration ratio pr/ pl should be equal to y Dl/Dr, 
instead of the required D^/Dr. The reason this argu- 
ment is wrong is that instead of a perfectly sharp jump 
in the concentration, it actually changes from pl to pR 
over a length scale on the order of the length of the par- 
ticle jump. Therefore the particle concentrations at the 
points where the jumps across the interface originate are 
not equal to the values far from the interface (pL and pr, 
respectively). Note that, while this is true even when r 
is continuously distributed, in the case r = ±1 there is 
an additional complication that affects the practical use 
of the BD algorithm: the particle concentration keeps 
varying wildly even far away from the interface, with the 
typical length scale of the oscillations again on the or- 
der of the particle jump size. As a consequence, the re- 
sults depend strongly on how the data are binned. For 
this reason, using r from a discrete set is best avoided in 
practical applications of the BD algorithm. 
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V. SIMULATION PROTOCOL 

The system is defined such that the left side extends 
from a wall at x = —d/2 to the interface at x = and 
the right side extends from the interface to a wall at d/2, 
where d is the total length of the system. Thus the in- 
terface is placed equidistant between the walls. For all 
simulations described in this paper, the particle is ini- 
tially placed at the interface. In the case of the MC II 
algorithm, when there is no site exactly at the interface, 
it starts at one of the sites immediately to the left and to 
the right of the interface with equal probabilities. Two 
types of simulations are performed. First, simulations 
with reflecting walls at x = —d/2 and x — d/2 are con- 
ducted. To obtain the stationary particle distributions, 
relatively long simulations must be performed and thus 
we focus on a single case where rjn/rji, = 4. Additional 
simulations with absorbing walls are performed to inves- 
tigate whether particles preferentially end up at one wall 
(herein referred to as an ultimate preferential direction) 
and to measure the conditional mean first passage time 
(MFPT) to each wall. In most of these simulations, the 
viscosity on the left-hand side, r\L, is held fixed and the 
viscosity on the right-hand side is varied from rjR = O.lrjL 
to 7]r = 10?7l; we also use tjr/tjl = 50 or even, in Ap- 
pendix A, r\R/r\L — 1000. 

For static quantities, such as the stationary particle 
distributions and the ultimate preferential direction, only 
the ratio of the viscosities across the interface influences 
the results — the actual values of the viscosities and the 
proportionality factor k in Eq. ([3]) do not matter. Indeed, 
for instance, looking at the expressions for the probabili- 
ties of the moves of the MC I algorithm [Eqs. (|5"8"]l. (1ST)1) . 
(|68|) . (|70|) . and flTTJ], if t]l and rjn are both changed by 
some factor and at the same time the adjustable constant 
?7q is changed by the same factor, then the probabilities 
of the moves remain the same. What does change is the 
time scale, since changing rj\ changes the time step tj, 
according to Eq. (|5D|) , which corresponds to speeding up 
or slowing down the dynamics. Since all approaches are 
expected to also give dynamically correct results such 
that the MFPTs and the decay of the transient behavior 
towards the stationary state should be consistent for any 
particular a value (with BD corresponding to a = and 
LD to a = 1), a comparison of the dynamical quanti- 
ties between the approaches is also meaningful. To make 
such comparisons, it is obviously necessary to know how 
the time steps of different algorithms are related. 

It is convenient to choose the time step of the MC I 
algorithm as the unit of time and express all other times 
and time steps in terms of this unit. That is, numerically, 
the time step of MC I will be chosen to be 



n = l. 



(86) 



Likewise, the mesh step of MC I is chosen as the unit 
of length (ai = 1). In these units, we use the system 
size d = 50, except for a single case in Appendix A. 
Thus there are 51 lattice sites at integer values of x from 



x = —25 to x = 25, and the interface at x = coincides 
with one of the sites as appropriate. The viscosity is 
expressed in units of rj^; thus, numerically 
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According to Eq. (|6T))) . in these units 



2kT 



1. 



In MC II, we choose arbitrarily aim = v2ai = v2 
in Eqs. (|T6"|) and ([77]). The lattice is placed in such a 
way that the interface is at distance a\ Y /2 from the site 
adjacent to it on the left and at distance a^/2 from the 
site on the right. For the leftmost site, —d/2 < x < 
—d/2 + a^, and for the rightmost site, d/2 — < x < 
d/2. It is convenient to express the viscosity in units 
of ryj 1 . However, for direct comparison with MC I the 
viscosity unit should be the same in both cases; thus we 
need r$ = rj\ (=1). Using Eqs. ([BDJ and ([751). 
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or, using the chosen values for auo and ai and keeping in 
mind that we have chosen the rest probability in MC II 
piis = 2/3 (as explained in Sec. HVBj) . 
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While choosing the values of the viscosities, we need 
to keep in mind that in MC I the viscosities cannot be 
smaller than t]q — 1. For the reflecting wall simulations, 
•qi, is set to 2, since, as we have mentioned, using a value 
too close to unity can introduce artifacts in the distribu- 
tion function. Correspondingly, ijr = 8 for r\R/r\h = 4. 
For the absorbing wall case, for the simulations where we 
vary r\R/r\h from 0.1 to 10, t]l is set to 10 so that r\R can 
be varied from 1 to 100. Using r\R ■ = 1 is not a prob- 
lem in this case, since we are not interested in the details 
of the particle distribution or the FPT distribution, but 
only the MFPTs, in which case all oscillations and other 
artifacts average out. For r\R/r\h = 50 and 1000, we use 
rjL = 2, as in the reflecting case, so rjR is respectively 100 
and 2000. 

For both MC algorithms, when the particle attempts 
to move left from the leftmost site or right from the right- 
most site, in the case of absorbing boundaries it disap- 
pears and the simulation is stopped, while in the case of 
reflecting boundaries the move is rejected. 

As for the MD simulations (BD and LD), we likewise 
use d = 50 and treat the particle as pointlike in the sense 
that it can approach the walls infinitely closely. The ran- 
dom terms in the equations of motion are normally dis- 
tributed. When attempting to cross a wall, the particle 
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disappears in the absorbing case and is "reflected" in the 
reflecting case. The latter means that, for example, when 
the particle crosses the right wall at d/2 and would end 
up at x = d/2 + A, it is placed at x = d/2 — A instead and 
in the LD case the sign of the particle velocity is changed; 
similarly for crossing the left wall. The equations of mo- 
tion that are solved numerically [Eqs. (fTCTj) and (j52")l . re- 
spectively] contain the diffusivity D. For comparison to 
MC, it needs to be expressed in terms of the viscosity 
given in the same units as in MC simulations, that is, in 
units of rj\. Using Eqs. ((SJ) and (|88p . in these units 
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For the BD time step St, it is reasonable to choose a 
value similar to the time resolution of the MC schemes. 
We have chosen St — 2 for simulations with reflecting 
boundaries and St = 20 for simulations with absorbing 
boundaries. While the latter may appear too large, we 
have checked that using a ten times smaller step pro- 
duces equivalent results. Ultimately, this is not surpris- 
ing: What matters is that the time step is much smaller 
than the typical time that it takes the particle to reach 
the boundary starting at the interface, and this condi- 
tion is satisfied. The situation with LD is somewhat 
more complicated. First of all, the initial velocity of the 
particle needs to be specified; it is chosen equal to ei- 
ther + y/WJ m or — ^JkT/m at random. Second, there 
is one other relevant combination of parameters, m/kT, 
and the associated time scale, mD/kT — m/2rjkT. It 
determines the inertia in the system, in particular, the 
correlation time of the particle velocity (or the time of 
relaxation to the terminal velocity). Since the intent is 
to model an overdamped system in which inertia is negli- 
gible, this time scale should be small, in any case smaller 
than the time it takes the particle to reach the wall from 
the interface. We chose m/kT = 4 for reflecting bound- 
aries and m/kT — 400 for absorbing boundaries. In the 
least favorable case, r\ = 1 for absorbing boundaries, the 
relaxation time is t r = m/2r\kT = 200. The typical par- 
ticle displacement over this time is \J2Dt r = \Jt r /r\ ~ 14 
which is less than a factor of 2 smaller than d/2 = 25. 
While some discrepancies between LD and the MC and 
analytic results are indeed found at the lowest viscosities 
studied, good agreement is obtained for the majority of 
rjR values. The time step in LD should be much smaller 
than the relaxation time; we choose 8t — 0.02 for reflect- 
ing boundaries and St — 0.2 for absorbing boundaries. 
The least favorable case is the largest r\ in our series for 
absorbing boundaries, rj = 100 (in the separate runs with 
7^=2000 we do not use LD, only MC I). The relaxation 
time in this case is t r — 2, so this condition is satisfied. 
Note that the time step for LD is much smaller than for 
BD and in fact, LD is the most time-consuming algorithm 
of all those considered in this paper. This is expected, 
since in LD the time step should be much smaller than 
the relaxation time, which should itself be small; we note, 
however, that the efficiency of the LD algorithm could be 



improved without sacrificing its accuracy by using differ- 
ent m/kT and St for different viscosity ratios. 



VI. RESULTS 

A. Reflecting Walls: Distributions 

In this section, we look at the particle concentration 
distributions for the case when the particles start at the 
interface. We first present the distributions for the Ito 
and isothermal calculi using the MC I and MC II algo- 
rithms, comparing them with the results obtained using 
BD (for Ito) and LD (for isothermal). We then compare 
the distributions for several different values of a using 
one of the MC algorithms (MC I). In all cases, we show 
the distributions at a time shortly after the start of the 
simulation and at a long time corresponding to the sta- 
tionary state. In addition, for the Ito and isothermal cal- 
culi we also include one intermediate time, such that the 
distribution on the low- viscosity side has already reached 
the quasistationary state, but that on the high-viscosity 
side still resembles the short-time distribution. For each 
curve, the total number of runs is 100000. In each run, a 
single particle is simulated and its position is determined 
at the times of interest. In MD simulations, the positions 
are binned with a unit bin size and the centers of the bins 
at half-integer x. The plotted values pi are proportional 
to the numbers of particles rii found at the given time 
at the site (or in the bin) i and inversely proportional to 
the distance between the sites (width of the bin). The 
distributions are normalized so the area under the curve 
is unity (thus, in effect, these are probability density dis- 
tributions). In particular, for MC I, BD and LD, since 
the distance between the sites (or the width of the bins) 
is unity, the normalization condition is 



Pi = !. 



where the sum runs over all sites (bins), and thus 



Pi 



For MC II, 
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where the first sum is over the sites in the left region and 
the second sum is over the sites in the right region, so 



il ( a ii Ej n j ) m the left region, 
if ( a nEj 71 j ) m tne ri g nt region. 



(96) 



We also plot the ratio of the numbers of particles on 
the right and left sides as a function of time. When 
calculating this ratio, particles exactly at the interface 
(which is possible in the MC I case) are ignored. 
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1. Ito calculus 

The probability distributions p as functions of position 
at three different times are shown in Fig. 03a)-(c), for all 
algorithms yielding Ito calculus (MC I, MC II and BD). 
In Fig. [3][a), the distribution shortly after starting the 
simulations, at time t — 100, is shown. By this time, only 
very few particles reach the walls, so the distribution is 
essentially the same as it would be in infinite space, with- 
out the walls. In this case, theoretically the distribution 
is expected to follow 
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This satisfies the diffusion equation P5|) in both regions, 
as well as the jump condition (|44|) and the flux matching 
condition (|4l))) . the latter because dp/dx is zero both at 
x = —0 and x = +0. Thus the parts of the distribution 
on the left side (L) and the right side (R) should both 
have a Gaussian form, but there should be a significant 
discontinuity in the distribution with a higher concen- 
tration of particles on the high viscosity side. These fea- 
tures are indeed reproduced in Fig. [3ja) . The disconti- 
nuity develops instantaneously (in a simulation, within 
a few steps from the beginning) and results from the 
trapping effect of high viscosity regions in the Ito for- 
mulation. Figure EJc) shows the distribution after a long 
time (t = 30000) thus corresponding to the stationary 
state. While the transient piecewise-Gaussian shape has 
disappeared, the discontinuity persists yielding a step- 
function form. Recall the Ito stationary state condi- 
tion [Eq. (|39l) ] which, using D ~ 1/rj, can be written 
as p(x)/rj(x) = const. Hence, as expected, the probabil- 
ity density on L, where the viscosity is low, is less than 
on R, where the viscosity is high, and the ratio is 



Pr/pl = Vr/vl- 
Given the normalization condition (p^ + pn)d/2 
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= 1, 
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(100) 



For the case presented in Fig. [3] where d = 50, tjl = 2, 
rjR = 8, we have p^ — 0.008, pr — 0.032. These values 
are indicated on Fig. (3Jc) where Pr/pl = 4 as expected 
for r\B.lr\L = 4. Note that at short times [as in Fig.^a)], 
the ratio of the probability densities immediately to the 
right of the interface and immediately to the left of the 
interface is also the same. In fact, this is true at all 
times, as follows from Eq. (T4"4"|) . The sharp drops in the 
distribution near the walls in Fig. [3Jc) are discretization 
and binning artifacts. 

In Fig. GDJd), the ratio of the number of particles on R 
to L, denoted Nr/Nl, is plotted as a function of time. 
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FIG. 3: (color online) Results of reflecting wall simulations 
corresponding to the Ito calculus, for r)L — 2, r/H = 8 and 
d — 50: particle concentration distributions, (a) at a short 
time, t — 100, (b) at an intermediate time, t = 1500, and 
(c) at a long time, t = 30000, as well as (d) the ratio of the 
number of particles on the right and left sides of the interface 
as a function of time, with the initial (t\) and final fa) times 
of the transition period marked. All particles start at the 
interface in the middle of the system. In all plots, the results 
obtained using both MC methods, as well as BD simulations, 
are shown. 



This ratio is equal to the ratio of the areas under the 
parts of the distribution curves on R and L. Note that 
Nr/Nl increases from a value around 2 and saturates at 
a value around 4. Hence, there are always more particles 
on the high viscosity side: Under Ito conditions, the high 
viscosity side acts as a trap. While the saturation value 
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agrees with the stationary state condition of Ito calculus, 
it is interesting that, beyond extremely short time be- 
havior which depends more on the simulation algorithm 
than the physics, the initial point for the increase during 
the transient portion is 2. Indeed, from Eq. valid at 
short times, the ratio should be \Jt\rIt\l = \/4 = 2. This 
is the ratio of the heights of the half-Gaussians (which is 
equal to Dl/Dr — r\Rl"r\h) multiplied by the ratio of 
their widths (which is equal to ^Dr/Dl = \Jr\hjriR)- 
At longer times, the particles begin to find the edge of the 
system and bounce back. Correspondingly, the system 
begins to move towards equilibrium where four times as 
many particles are found on the high viscosity side. The 
transition starts when the particles on the lower- viscosity 
side reach the boundary, which happens around time t\ 
such that 2Z? max ti = (d/2) 2 , thus ti = rf 2 /(8D max ), or, 
using Eq. O, 



d 2 r] n 



(101) 



where Z? max is the higher of the two values of the diffu- 
sivity and rj m i n is the corresponding (lower of the two) 
value of the viscosity. The transition ends when the par- 
ticles on the higher- viscosity site also reach the boundary, 
which occurs around 

d //max 



(102) 

For our values of the parameters (d = 50, r/min = t/l = 2, 
?7max = VR = 8), ti = 1250 and t2=5000. These two 
times are marked on Fig. G^d). The time at which the 
distributions in Fig. Eta) are obtained is much smaller 
than t±, thus corresponding to the situation before the 
transition starts, and the value of the ratio Nr/Nt, at 
that time is indeed close to two, according to Fig. EJd). 
On the other hand, the time at which the distributions 
in Fig. [3tc) are obtained is much larger than £2 and thus 
corresponds to the stationary situation after the transi- 
tion has ended; the value of Nr/Nt, is close to four. 

In Fig- EJb), we show the distribution at time t = 1500, 
which is between t\ and t^. By time t such that 



ti < t < < 2 , 



(103) 



particles on the low-viscosity side would have traversed 
that side many times, so the distribution on that side 
should be nearly flat resembling the stationary distribu- 
tion, but on the other hand, very few particles would have 
reached the wall on the high- viscosity side, so the distri- 
bution on that side should still be close to the short-time 
Gaussian distribution. Even though in our case the ratio 
t2/ti — r\Rlr\h is only four, so, strictly speaking, sat- 
isfying Eq. (| 1031) is impossible, Fig. Efb) can still serve 
as a reasonable illustration of that regime. Note that 
the distribution on the low- viscosity side is not truly sta- 
tionary, but rather quasi-stationary, since there is a net 
outflow of particles to the high-viscosity side, as indi- 
cated by the gradual increase of Nr/Nl between ti and 



t2 [see Fig. Eld)]. This outflow is slow enough for the 
distribution to remain nearly flat, but the height of the 
distribution gradually decreases. In essence, the height of 
the plateau on the low- viscosity side adiabatically follows 
the height of the peak on the high-viscosity side (which 
decreases as this peak broadens), so that the boundary 
condition ([4"4"]l is satisfied at all times. On the other hand, 
this outflow is large enough that the slope of the distribu- 
tion function is significantly nonzero immediately to the 
right of the interface in Fig. Etc). However, the viscosity 
ratio is rather small here and also t is close enough to 
t2 that particles have already started to reach the right 
wall (note a nonzero particle concentration next to it). 
As we argue in Appendix A, for a very high viscosity ra- 
tio and when the condition (|103[) is truly satisfied, the 
derivative at x = +0 is negligible and the distribution on 
the high- viscosity side is a Gaussian. 

Taking into account the general character of the dis- 
tribution at intermediate times between ti and t2, we 
can also obtain approximately the time dependence of 
the ratio Nr/Nl in this time interval. The num- 
ber of particles on R, Nr, is roughly proportional to 
p{+0, t) times the width of the distribution on that side, 
which is oc \ft. The number of particles on L, Nl, is 
proportional to p{— 0,t) times the length of that part 
of the system, d/2, which is time-independent. Thus 
Nr/Nl oc [p(+0, t)/p(— 0, t)]y/i, and, given that, accord- 
ing to Eq. flUJ) , the ratio in the square brackets is time- 
independent, 

^ cx Vt. (104) 

This agreement between the theory at a = and the 
Brownian dynamics results confirms that the BD sim- 
ulations correspond to Ito calculus. Further, note the 
overall good agreement between the BD simulations and 
both MC approaches for all the results shown in Fig. El 
Considering factors such as discretization, binning, and 
noise, the fact that the differences are small between the 
data sets is a convincing validation of the two MC ap- 
proaches we have developed to simulate Ito calculus. 



2. Isothermal calculus 

The equivalent plots for all algorithms yielding the 
isothermal calculus (MC I, MC II and LD) are shown 
in Fig. 01 Note that the distributions are now continu- 
ous across the interface, in agreement with Eq. (|44"1) for 
a = 1. For the short time behavior shown in Fig. Ufa), 
where, again, t = 100, the Gaussian-like distribution on 
the left side is wider than that on the right side, as in the 
Ito case. Theoretically, the distribution at short times is 
given by 
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(105) 
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which is the same as Eq. (p?T| for the Ito case, except for 
the prefactors that ensure the continuity across the inter- 
face. At long times, in the stationary state, there is an 
equal probability to find a particle anywhere in the sys- 
tem, as is seen in Fig. @|c) (t = 30000). This is consistent 
with the isothermal stationary state condition [Eq. (j42|) ] . 
The stationary value of the density is obviously given by 



p=- = 0.02. 

a 



(106) 



The time dependence of the ratio of the number of par- 
ticles on R over L is shown in Fig. HJd). Here, the tran- 
sient behavior begins from the ratio of approximately 0.5. 
That is, given the bias at the interface to end up on the 
lower viscosity side, very rapidly there is a greater num- 
ber of particles on L, and, contrary to the Ito results 
where there are twice as many particles on the high vis- 
cosity side near the beginning, for the isothermal results 
there are half as many particles on the high viscosity side. 
This ratio can again be obtained by considering the peak 
height and width. Here, the peak height is the same on 
both sides while the ratio of the widths is again \J^]hh\k- 
Correspondingly, the ratio is Nr/Nl = y/rjE/ \/Wr = l/2< 
As the system moves towards equilibrium, Nr/Nl again 
climbs and plateaus, but now it saturates at 1.0 in accord 
with the isothermal stationary state condition. The tran- 
sition between the initial and the final values again starts 
around t\ given by Eq. (|101j) and ends around t% given by 
Eq. (I102p . At times between t\ and t-z , the distribution 
is nearly flat (quasi-stationary) on the low-viscosity side 
and Gaussian-like on the high- viscosity side, again, as in 
the Ito case [see Fig. [4f b) for t = 1500]. However, unlike 
in the Ito case, we show in Appendix A that even for a 
large viscosity ratio and t% <C t <C t%, the slope at +0 is 
significantly nonzero, and, in fact, the interface lies at the 
inflection point of the Gaussian-like curve rather than at 
its maximum. Despite this, the crude arguments leading 
to Eq. (|104p still apply, so we expect that time depen- 
dence to still be approximately valid here. We compare 
the time dependences of Nr/Nl for different calculi in 
the next section. 

In all four plots, again, there is good agreement be- 
tween all three algorithms. 
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3. Comparison between calculi 

In Fig. [51 we compare the results for the particle dis- 
tributions and the ratio Nr/Nl for several different cal- 
culi (values of a), including the already considered Ito 
(a = 0) and isothermal (a = 1) limit cases, using the 
MC I algorithm. Essentially identical results are obtained 
using MC II (not shown). Figure [5j(a) shows the distribu- 
tions at t = 100 and Fig. [5jb) presents the distributions 
at t = 30000. At short times [Fig. [SJ^a)], the distribution 



FIG. 4: (color online) Results of reflecting wall simulations 
corresponding to the isothermal calculus, for r\L — 2, t}r — 8 
and d — 50: particle concentration distributions, (a) at a 
short time, t — 100, (b) at an intermediate time, t = 1500, 
and (c) at a long time, t = 30000, as well as (d) the ratio 
of the number of particles on the right and left sides of the 
interface as a function of time, with the initial (ti) and final 
(£2) times of the transition period marked. All particles start 
at the interface in the middle of the system. In all plots, 
the results obtained using both MC methods, as well as LD 
simulations, are shown. 



18 



is always piecewise-Gaussian with 

D 

p(x,t) 



D a L - 1/2 +D«- 1/2 V^i 

DT 1 1 



exp 



4D L t y ' 



D" 



2 +D a ~ 1/2 v^* 



exp 



\ 4D R t J 



x < 0, 
x > 0. 
(107) 

There is a jump at the interface in all cases, except a = 1, 
in agreement with Eq. (|44p . At long times [Fig. [SJb)] , the 
distribution is piecewise-uniform, with 



PL 



PR 



(108) 
(109) 



which follows from Eq. (|37j) for the stationary distribu- 
tion. Again, there is a jump for a ^ 1. In Fig. [Sfc), 
the ratio Nr/Nl is plotted. The short-time value of this 
ratio is 



Nr 
N L 



l/2-a 



(110) 



which is the product of the ratio of the peak widths 
[(vr/vl) -1 ^ 2 ) and the ratio of the peak heights 
[(Vr/Vl) 1 ^ 01 }- The long-time (stationary-state) value of 
N R /N L is 



Nr 



Vr 

VL 



(111) 



For < a < 1/2, Nr/Nl is always larger than 1, so there 
are always more particles on the higher-viscosity R side 
than on L. For a = 1/2, the ratio is one initially, but then 
grows above unity. For a = 1, there are always more par- 
ticles on the lower-viscosity site, but the ratio approaches 
unity from below as t — > oo. Finally, for 1/2 < a < 1, 
there are more particles on the low- viscosity side initially, 
but the ratio crosses unity at a finite time and eventu- 
ally there are more particles on the higher-viscosity side. 
The transition between the short-time and the long-time 
values of Nr/Nl is roughly between ti and t 2 given by 
Eqs. (jlOip and (|102p . and the time dependence should 
still be roughly given by Eq. (|104l) . 

Note that at both small and large times, if the parti- 
cle distribution on L is multiplied by (t/r/vl) 1 ^ 01 , then 
the resulting function is the same for all a: This rescal- 
ing eliminates the jump at zero, and the shapes of the 
distributions in a given region and at any given time are 
calculus- independent, both at short times (when they are 
Gaussian with the calculus-independent width) and at 
long times (when they are flat). A reasonable question 
then is if this is also true at intermediate times and thus 
at all times. The answer, perhaps somewhat surprisingly, 
is no. This can be seen by comparing the shapes of the 
intermediate-time distributions for the Ito [Kig. |3jb)] and 
isothermal [Fig.HJb)] cases that are indeed somewhat dif- 
ferent, and, as mentioned, the difference is significant in 
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FIG. 5: (color online) Results of reflecting wall simulations 
using the MC I algorithm for different calculi, with the values 
of a specified in the legend, for vl — 2, nn — 8 and d = 50: 
particle concentration distributions, (a) at a short time, t = 
100, and (b) at a long time, t = 30000; (c) and (d) the ratio of 
the number of particles on the R and L sides of the interface 
[rescaled in (d) as indicated] as a function of time, with the 
initial (£i) and final (£2) times of the transition period marked. 
The rescaling factor in (d) is chosen so that the values coincide 
at both small and large times, yet in the intermediate region 
small differences are observed. The straight line indicates the 
predicted ~ t 1 ' 2 behavior. All particles start at the interface 
in the middle of the system. 
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the high- viscosity-ratio limit (see Appendix A) . This dif- 
ference is also easy to explain: If we take, for example, 
the solution for the Ito case and multiply its left half 
by rin/riL, the resulting function will satisfy the diffusion 
equation (|45[) and the jump condition (|44l) for a = 1, but 
not the flux continuity condition (|46p , since only the left- 
hand side of the latter will be multiplied by rjn/rjL, but 
not the right-hand side. The cases t — > and t — > oo are 
exceptions, because in these limits dp/dx is zero at both 
x = — and x = +0. However, for intermediate times 
this cannot be the case, since there must be some net flux 
between the two sides, or otherwise Nr/Nl would remain 
constant. As a consequence, while the curves in Fig.(5Jc) 
for different a do look similar, and, in particular, the 
transition between the low-time and the high-time values 
of Nr/Nl occurs over roughly the same time interval and 
according to approximately the same ^ t 1 / 2 dependence 
in all cases, these curves are not perfect rescaled copies 
of each other: Multiplying these functions by (r)n/riL) a 
so they match at t — ¥ and t — ¥ oo does not make them 
coincide everywhere, as illustrated in Fig. [5jd), where 
small differences between the rescaled curves around t2 
are apparent. That plot also compares the curves with 
the t x l 2 dependence predicted for large viscosity ratios; 
the agreement is reasonable given that in our case the 
ratio rjn/ijL = 4 is not large. 



B. Absorbing Walls: Preferential Direction 

We now describe the results of simulations with ab- 
sorbing walls, that is, each single-particle simulation 
run is stopped once the particle reaches one of the two 
walls. For the Ito and isothermal calculi, we first look at 
whether the particles prefer to end up on one of the walls 
(i.e., whether there is an ultimate preferential direction) 
and find the conditional MFPTs, separately for parti- 
cles reaching each wall. The ultimate preference is then 
compared to the preference of the particle to spend more 
time on one or the other side, which is analyzed in two 
different ways. Unlike the case of reflecting boundaries, 
where a single set of viscosities was considered, here we 
vary the ratio of the viscosities, with the viscosity on L 
kept constant. For each set of viscosities, we carry out 
10000 single-particle runs, again starting at the interface, 
as described in Sec. |V] For the Ito calculus we also ob- 
tain the full conditional FPT distributions for a single 
representative viscosity ratio using 10 6 runs. 



1. Ito calculus 

Results for the absorbing wall simulations in the case 
of Ito calculus are shown in Fig. [BJ Figure [SJa) displays 
the fraction of events absorbed on R, itr (dashed lines) 
and L, ttl (solid lines) for all three algorithms. For each 
algorithm at any viscosity ratio, the number absorbed at 
L equals the number absorbed at R within error. Thus, 



under Ito conditions, there is no ultimate preferential di- 
rection in the system. This result may seem paradoxical 
at first, given that for large rjn/r/L there are much fewer 
particles on L than on R at any given moment of time, 
which we have seen for reflecting boundaries [Fig. [31(d)] 
and which turns out to be true for absorbing boundaries 
as well. The primary reason for the equality in Fig. [S] 
(a) is that particles diffuse faster on L, so even a smaller 
concentration creates a significant flux towards the left 
wall. 
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FIG. 6: (color online) Results of absorbing wall simulations 
corresponding to the Ito calculus, for vl = 10, d = 50, and 
varying vr: (a) the percentage of particles absorbed by L 
(solid lines) and R (dashed lines) walls; (b) the MFPT under 
the condition that the particle is absorbed by the left wall 
(solid lines) and by the right wall (dashed lines) . All particles 
start at the interface in the middle of the system. In both 
plots, the results obtained using both MC methods, as well 
as BD simulations, are shown. Black lines in (b) are linear 
fits to the averages of all three methods. 

The conditional MFPTs to each wall (which we de- 
note Tl and tr) are plotted in Fig. EJ^b). As would be 
expected, the time to the lower viscosity side is always 
less than the time to the high viscosity side with an in- 
tersection corresponding to equal times at ijn/r/L = 1. 
A less obvious fact is that the two times are always of 
the same order of magnitude even when the viscosities 
are very different. In fact, both dependences are lin- 
ear as demonstrated by straight-line fits (black lines) to 
the averages of all three algorithms, and the slope for 
the R data is, within the uncertainty, two (2.04 ± 0.07) 
times that of the slope for the L data, so the ratio of 
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the two times approaches two as rjn/rjL — > oo. This may 
seem surprising, given that the time t\ [Eq. f|101[) ] around 
which a significant (not exponentially small) number of 
particles first reaches the left wall is much smaller in this 
limit than the corresponding time t^ [Eq. (|102|) ] for the 
right wall. The explanation is that, while it is indeed true 
that the particle flow to the wall on the low- viscosity side 
starts at time ~ t\, it continues until time ~ £2, as the 
particles keep "leaking" from the high-viscosity side un- 
til those that have not been absorbed by the wall on the 
low-viscosity side are absorbed by the other wall. Even if 
the means are of the same order of magnitude, the con- 
ditional first-passage time distributions PL,n(t) are very 
different when the viscosities are very different. For the 
wall on the low- viscosity side, the flux reaches its max- 
imum very fast, by time ~ t\ and then decreases very 
gradually (as a power law) until time ~ £2, after which 
it decays exponentially. But for the wall on the high- 
viscosity side, the flux is essentially zero until time ~ £2 
and then rises and decays relatively quickly. This behav- 
ior of the distributions is illustrated in Fig.[7]obtained for 
a large viscosity ratio {t\l = 2, t]r = 100, so r]n/r]L = 50) 
using the MC I method with the waiting time selected 
from a distribution as given by Eq. (|72")l . The best fit to 
the power-law region has the exponent pa —0.8; we show 
in Appendix A that in the limit of a very large viscosity 
ratio this exponent is —0.5. 




FIG. 7: (color online) The conditional FPT distributions for 
the left and right walls for r/L — 2, r)R — 100, d — 50 obtained 
using the MC I algorithm. Times ti and ti given by Eqs. (|10ip 
and (|102p . respectively, are marked. 



To derive the properties revealed in Fig. [6j consider the 
variant of the MC I algorithm with all steps successful 
(no staying put) and the time increment at each step t w 
given by Eq. ([73} (see Sec. IIVA2|) . Recall that in the 
Ito formulation, there is no bias at the interface [from 
Eqs. dTJJl and ([75]). pf_ = p[° + = 1/2], or anywhere else 
in the system. Hence, the trajectory of a particle is just 
an unbiased random walk, except for the fact that the 
time increment at each step changes depending on what 
side of the interface the particle is at. This exception, of 
course, does not change the fact that this random walk is 



equally likely to first reach the left or the right wall, since 
it starts at equal distances from them, which explains the 
result in Fig.^a). To get insight into the results for the 
MFPTs in Fig. \$[b), we use the fact that if a particle 
is released equidistant between two absorbing walls that 
are at distance d from each other and this particle does 
a random walk with step size unity eventually reaching 
the right-hand wall, then on average it spends 
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steps on the left-hand side (see Appendix B). Then 
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r[i+2(m/vL)], (H4) 



■l = ti w n same + t£,n opp = ^^[2 + {r) R /r) L )}, (115) 



where t\ w and t^ w are the time steps [given by Eq. ([73]) 
on L and R, respectively, or, using the conventions ([86 
and (15711 . 
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For til — const, the dependence of both times on Jjn/rii 
is indeed linear, and the ratio of the slopes is two, in 
agreement with the data in Fig.[6](b). In our case [d = 50, 
VL = 10), 



tr ps 2083 + 4167(7?*/^), 
t l « 4167 + 2083(7^/7^). 



(118) 
(119) 



This should be compared with the fits to the data (black 
lines in Fig. [B]), which are 



Jit 



Jit. 



(2.35 + 0.03) x 10 3 
+ (4.26 + 0.09) x 10 3 (W?7l), 
(4.58 + 0.04) x 10 3 

+(2.09 ± 0.06) x 10 3 (r te /77 L ). (121) 



(120) 



The slopes are indeed in agreement; the small disagree- 
ment in the intercepts (about 10%) is likely due to dis- 
cretization effects. In the limit rjR/r/L — > 00, Eqs. (|116[) 
and (|117[) become 
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where £2 is given by Eq. (|102p . This shows that in this 
limit both times become independent of the lower vis- 
cosity and in particular, both are ~ ti. We discuss the 
relation between this result and the FPT distribution in 
Appendix A. 



2. Isothermal calculus 

The same data for the case of isothermal calculus are 
shown in Fig. [8j There are significant quantitative and 
qualitative differences compared to the case of Ito calcu- 
lus. First, examining Fig. IFJa), there is a clear ultimate 
preferential direction in the system. When t\rIt]l = 0.1, 
approximately 90% of the particles are absorbed at the 
R wall. Conversely, when ijr/tjl — 10, 90% of the parti- 
cles are absorbed at the L wall. Plotting the data with 
a log scale for the viscosity ratio [inset to Fig. |8ja)] 
demonstrates that this symmetry holds for all simulated 
viscosity ratios: 



(124) 



where ttl.r{ z ) are the probabilities of absorption at the 
left (right) wall when the ratio of viscosities on R and on 
L is z. This symmetry is a necessary outcome considering 
that in the algorithms, it is only the ratio of viscosities, 
and not the absolute value, which influences the dynam- 
ics (beyond simply making all events take longer). The 
ultimate preferential direction under isothermal condi- 
tions arises from the bias at the interface favoring the 
low viscosity side. Since this bias increases with a larger 
disparity in viscosities, the strength of the preference in 
direction also grows. It must, however, saturate at 100% 
resulting in the behavior shown in Fig. [8ja). 

The MFPTs, shown in Fig. [5Jb), are also quite differ- 
ent from those corresponding to the Ito condition. Again, 
as expected, the time to the lower viscosity side is lower 
than the time to the high viscosity side with an inter- 
section at rjn/riL = 1. However, the dependence across 
the entire viscosity range is clearly not linear. In fact, 
considering data for which rjn/rjL > 1, while tr does in- 
crease linearly with increasing ?/.r/?7z, (albeit only when 
this ratio is large), 77, appears to saturate. 

To understand the results in Fig. [5J we turn again to 
the variant of the MC I algorithm with all steps success- 
ful. For generality, we consider an arbitrary a and then 
specialize to a = 1. For a/0, the segments of the parti- 
cle trajectory between visits to the interface site are still 
random walks, but when the particle is at the interface, 
it chooses to continue on L or on R with the probabili- 
ties given by Eqs. (|74l) and (|75|) . respectively. Each such 
trajectory segment is contained entirely either within the 
left half or within the right half and for brevity, we refer 
to these groups of segments as left and right segments. 
Since the probabilities to reach the respective wall during 
the left and right segments are equal, the total probabil- 
ity to reach the left wall, ttl, is equal to the fraction of 
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FIG. 8: (color online) Results of absorbing wall simulations 
corresponding to the isothermal calculus, for t\l = 10, d = 50, 
and varying t/r: (a) the fraction of particles ttl,r absorbed 
by L (solid lines) and R (dashed lines) walls; (b) the MFPT 
under the condition that the particle is absorbed by the left 
wall (solid lines) and by the right wall (dashed lines). Black 
lines in both plots are theoretical dependences derived in the 
text. All particles start at the interface in the middle of the 
system. In both plots, the results obtained using both MC 
methods, as well as LD simulations, are shown. 



left segments, which in its turn is equal to the probabil- 
ity to start a left segment while at the interface, which is 
p'i_. Correspondingly, the probability to reach the right 
wall, 7Tr, is p?,. So 
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7Ti? = 
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(125) 
(126) 



Note that — > 1 for i]r/t]l — > 00 whenever a > 0, 
whereas for a = 0, ttl = 1/2, as we have seen in the 
previous section. In the isothermal case (a = 1) 
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These functions are plotted in Fig. [5Ja), and agreement 
with the simulation data is excellent. Equations (|127[) 
and (|128|) satisfy the symmetry condition (|124j) . At 
Vr/Vl = 0.1, ttr = 10/11 ps 91%, while at VR /r] L = 10, 



22 



til = 10/11 f=a 91%, again, in agreement with the simu- 
lations. 

To calculate the MFPTs, we use the result obtained 
in Appendix B that for a random walk that is unbiased, 
except at the interface (where it chooses to go left with 
probability or right with probability the mean 
number of steps that the walk absorbed by the right wall 
spends on R is 



n RR = — (l + 2p&), 
while on L it spends on average 
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for a — 0, 

for < a < 1, (138) 
for a = 1. 



For = const, this always gives a linear dependence 
on (vr/vl) for tr, regardless of a, and, unsurprisingly, 
tr ~ t2, again for any a. For tl, the dependence is 
oc (vr/vl) 1 ^ 01 , which is linear for a = 0, reproducing the 
result of the previous section, sublinear for < a < 1, 
and a constant for a = 1, in agreement with Fig. [8f b) . 
Note also that tl is now of order t1t\~ a . This is deter- 
mined by the power-law region between t\ and t^ in the 
FPT distribution that still exists for a > 0, as discussed 
in Appendix A. 



steps. Then the MFPT for the walk ending on the right 
wall is 



T R — t\w n RL + tfw n RR 
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where we have used Eq. ([73")) for tj w and tf w and conven- 
tions (1551) and (|87p. For walks ending on the left wall, 
analogously there arc on average 



steps on L and 
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steps on R, and then 



3. Preferences during the diffusion process 

It is interesting to compare the preferences of the par- 
ticles to be absorbed by a particular wall (the ultimate 
preferential direction at the end of the trajectory) to their 
preference for a particular side at earlier times. The prob- 
lem can be formulated in different ways. 

First, we can take a sample of particles, calculate the 
total combined time spent by them on a particular side, 
and divide by the total combined duration of all trajec- 
tories. Obviously, this is the same as the average time 
that a trajectory spends on a particular side divided by 
the total average duration of the trajectory. Using again 
the MC I algorithm with all moves successful, we note 
that the average number of steps in each segment is the 
same for left and right segments. Denoting the average 
total number of steps n and given that the fraction of left 
segments is pf_ , as discussed in the previous section, the 
average time a trajectory spends on L is 
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In particular, in the isothermal case (a = 1), 
d 2 rj L (r]R/r] L ) 2 + 5{n R /r] L ) 
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Equations (flgS]) and (ITM| with r\ L = 10 and d = 50 
are plotted in Fig. HJb). The agreement with the sim- 
ulations is very good; discrepancies can be attributed 
to discretization effects as well as inertial effects in the 
LD case. Equations (|116|) and (|117|) are recovered from 
Eqs. (|131|) and (|134j) for a — 0. When t]r/i]l is large, 
Eqs. ([TBTj) and ATM]) give 




for a = 0, 
for < a < 1, 
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similarly the average time on R is 
t R = np!? + t? w , 
and then the fractions of the time on L and R are 
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These functions for a — 0, 0.5, and 1 are plotted in 
Fig. [HJa) together with the results of MC I simulations 
in these cases, and agreement is excellent. Note that 
these probabilities are different from the probabilities to 
reach a particular wall, ttl,r, given by Eqs. (|125p and 
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(|126[) . In particular, while for Ito ttl — ttr, /j, < f R . 
when tjl < t]r. Conversely, for isothermal /j, = f R , but 
TTi > ti\r for ?7l < rjR. In fact, 7Tl for Ito equals f R for 
isothermal, and vice versa, or, in general, 

fn(VR/VL,a) = ir L (n R /n L , 1 - a), (143) 
fL{vn/VL,a) = n R (r) R /r)L,l- a). (144) 

At the same time, /l,r coincide with the long-time prob- 
abilities for particles to be found on a particular side in 
the problem with reflecting boundaries [cf. Eq. (Ill ip ] - 
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9: (color online) For the case of absorbing walls with 
r)L = 10, d — 50, and varying rjn, for Ito, Stratonovich, and 
isothermal calculi, fractions of time spent by trajectories on L 
and R defined in two different ways, as described in the text: 
(a) }l,r\ (b) 4>l,r- Black lines are analytic dependences de- 
scribed in the text: Eqs. (|141[) and (|142[) in (a) and Eqs. (|145j) 
and (|146[) in (b). In all cases, particles start at the interface 
in the middle of the system. 

Second, we can also determine the fraction of the time 
spent on a particular side for each trajectory in the sam- 
ple and then average these fractions. We denote these av- 
eraged fractions (\>l,r- This definition gives more weight 
to short trajectories than the previous one, since a short 
trajectory contributes little to both the total duration 
of trajectories and to the duration on a particular side, 
the quantities that enter the definitions for Jl,r, but it 
contributes as much as long trajectories to 4>l,r- The 
results for 4>l,r. obtained in MC I simulations for a = 0, 



1/2, and 1 are shown in Fig. IHJb). These results can be 
fitted using the approximate expressions, 
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If the equalities (| 145|) and (|146j) were exact, then, (1) 
4>l for Ito would be equal to (f> R for isothermal, and (2) 
(f>L and (f>R for Stratonovich (a = 1/2) would be exactly 
1/2 for all r\nl r \L- However, we have checked using sev- 
eral different lattice sizes and extrapolating to the infinite 
size (results not shown) that these relations are indeed 
approximate in the continuum limit and both statements 
(1) and (2) are approximate as well. 

Comparing Eqs. (fT25|) . (fT26j) . (fTITj) . (fT42l . (jT45]) . and 
(|146|) . we note that the average fractions of time spent on 
the low- viscosity side calculated in both ways are always 
smaller than the probability to be absorbed by the wall 
on the low- viscosity side, but the second way of calculat- 
ing the average gives a result intermediate between the 
other two. That is, for r\ R jr\L > 1, 



II <4>L < 7T L - 



(147) 



The differences between these quantities are particularly 
striking when the viscosity ratio is large, in which case 
both "<" signs in Eq. (|147l) turn into so there are 

ranges of a where one or two of the quantities /z, and (f>L 
approach zero while itl approaches unity. These relations 
are easy to explain by recalling that many particles that 
are eventually absorbed by the wall on the low- viscosity 
side linger until time 3> t\ (up to ~ £2) necessarily spend- 
ing much of that time on the high- viscosity side, which 
explains why < til- In addition, such long tra- 

jectories have a relatively larger weight when calculating 
/l than when calculating 4>l, which explains /j, < 4>l- 
Note that nearly all trajectories reaching the wall on the 
high-viscosity side spend most of the time on that side. 



VII. DISCUSSION 

In this paper, we have developed two dynamical Monte 
Carlo algorithms for simulating the dynamics of a par- 
ticle in the presence of sharp viscosity (or, in general, 
diffusivity) changes applicable to cases where different in- 
terpretations of the stochastic term in the corresponding 
overdamped Langevin equation (or calculi, such as Ito, 
Stratonovich, and isothermal) are appropriate. The va- 
lidity of these algorithms was demonstrated by the good 
agreement between the approaches and with Brownian 
dynamics for the Ito calculus and Langevin dynamics for 
the isothermal calculus, as well as with theoretical predic- 
tions. We have argued that MC approaches, such as those 
developed here, can be particularly advantageous when 
a sharp interface is present in the system, since solving 
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the overdamped Langevin equation directly is problem- 
atic in this case, except for the Ito calculus (correspond- 
ing to BD), while the LD approach is only applicable to 
the isothermal case and generally more computationally 
intensive. 

While our simulations of particles on a ID interval 
with the interface in the middle were intended primar- 
ily to validate our approaches, the system itself turned 
out to be rather interesting, with sometimes non-obvious 
properties that depend significantly on the choice of the 
calculus characterized by the parameter a that can vary 
between zero (Ito) and one (isothermal) . The differences 
between these extreme cases are summarized in Table U 
for a number of features of the system. Not only do the 
features vary greatly between calculi, but there are un- 
expected symmetries between certain features (indicated 
by arrows in the table) . Clearly, when modeling a physi- 
cal system, the choice of calculus in deciding how to treat 
the interface must be done with great care. 

The goal for the future is extending the approaches pre- 
sented here to more complex and physically interesting 
situations, such as two- and three-dimensional systems 
and biased diffusion under the action of a force. We cur- 
rently work on applying these algorithms to studying the 
translocation of a polymer through a hole in a membrane, 
with fluids of different viscosities on the two sides of the 
membrane. Of course, extension to the case of several 
interfaces is straightforward and a gradual change in vis- 
cosity can most likely be modeled as many interfaces one 
lattice site apart. Besides this, even in the simplest case 
considered here, the algorithms can be further improved: 
There are several free parameters that we have chosen 
arbitrarily, but that can be optimized to speed up the 
convergence to continuum, so even coarse discretization 
produces accurate results. 

It is useful to note that the physical scenario examined 
in this paper actually belongs to a more general class 
of problems. Of course, the discrimination between dif- 
ferent calculi arises for any stochastically driven process 
that can be written as an overdamped Langevin equa- 
tion in which the prefactor of the noise term is state- 
dependent [1, H, El HI [H [H. Hence, while we have 
focused on diffusion in an inhomogenous medium, simi- 
lar concepts can be ap plied to many problems including 
financial modeling [5J, [55ll. NMR spectroscopy [13] , pop- 
ulation growth models [5q| . and climate modeling [58l — 
lotij . Another interesting example is the work by Smythe 
et al. [61( examining noise-induced phase transitions 
where it was found that digital simulations produced re- 
sults in agreement with the Ito formulation while appar- 
ently equivalent analog simulations produced results that 
agreed with the Stratonovich formulation. Given these 
diverse scenarios, the simulation approaches introduced 
herein may be useful outside of the problem for which 
they were originally intended. 
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Appendix A: Particle distributions and 
first-passage-time distributions at intermediate 
times in the high-viscosity-ratio limit 

To explain the power-law dependence of the FPT dis- 
tribution in the intermediate-time regime between ~ t\ 
[Eq. ([TOll ] and - t 2 [Eq. (|T02l) ]. as observed in Fig. 
consider the particle distribution function p(x, t) in this 
regime for the problem with absorbing boundaries. Sim- 
ilar to the case of reflecting boundaries discussed in 
Sec. I VI Al p(x,t) consists of a quasistationary distribu- 
tion in the low- viscosity region and a "Gaussian-like" dis- 
tribution in the high- viscosity region. (The exact mean- 
ing of "Gaussian-like" will be revealed later.) The differ- 
ence is that because of the absorbing boundary conditions 
[p(x,t) = at the walls], the quasistationary part is no 
longer a constant, but instead a linear function, 



p(x,t) « A(t)(l + 2x/d), x < 0, 



(A.I) 



where A(t) = p(—0,t). Here and in what follows we 
assume that the left half of the system is the low- viscosity 
region (rjL <C r]n). Equation (|A.ip provides a relation 
between the value of p(x, t) at —0 and its derivative: 



dp(x,t) 



dx 



ip(-o,t). 



(A.2) 



The slope, [dp(x,t)/dx]\ x= _ d/2 w [dp(x,t)/dx]\ x= _ , is 
related to the flux to the left wall and thus to the condi- 
tional FPT distribution for the left wall, 



Dl dp(x,t) 

7Ti 



dx 



2D L 

7T L d 



A(t). (A.3) 



From Eq. (|A.2I) . using the jump condition (|4~4")) and 
the flux continuity condition (|46[) . we obtain a similar 
relation on the high- viscosity side of the interface: 



dp(x,t) 



dx 





J- 




x=+a 


d 





P (+0,t). 



(A.4) 



For the diffusion problem in the high-viscosity region, 
this serves as a boundary condition (BC). BCs of this 
type, where the derivative is proportional to the value 
of the function, are known as Robin boundary conditions 
[62| . At times t <C £2 the right boundary can be ignored, 
and we can treat this problem as being defined on a semi- 
infinite interval, 



dp(x,t) d 2 p{x,t) 
= L) R - 



dx 2 



x > 0, 



(A.5) 
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Feature 



Stationary concentration 
distribution 

Ultimate preferential direction 

Transient preference averaged 
per total time 

Transient preference averaged 
per trajectory 

Mean first passage times 



Distribution of first passage times 
at large viscosity ratios 

Relative number of particles 
on either side 



Ito: a=0 



Higher concentration on the more 
viscous side 



Time to either wall increases with 
increasing viscosity on one side 

Times to each wall on the same 
order of magnitude 

P ~ t at intermediate times 



short time: N = (t]j/1 l ) 
long time: N E /N L = 



Isothermal: a=l 



Equal concentration everywhere 



Time to the low viscosity side 
saturates, time to high viscosity 
side continues to increase 

Times to either wall can be very 
different 

P ~ t 13 at intermediate times 



short time: N^/IS^ = Oi^/tIl)" 1 
long time: N R /N L = 1 



No preference 



High viscosity side favored 



X 



High viscosity side favored 



More particles absorbed at the 
wall on the less viscous side 

k No preference 



►Low viscosity side favored 



TABLE I: Results for selected features of the system for the two extreme a cases: a = (Ito calculus) and a = 1 (isothermal 
calculus). The arrows indicate effects of equal magnitude (or approximately so indicated by ~). 



with the boundary condition (|A.4|) . 

The boundary-value problem given by Eqs. (IA.5[) and 
(|A.4|) can be solved exactly in quadratures [63J. We will 
opt instead for a more intuitive approach, and from now 
on, we will omit all numerical constants of order unity 
from our considerations. Since the BC (|A.4j) is interme- 
diate between the absorbing BC p(+Q,t) — and the 
reflecting BC [dp(x,t)/dx]\ x=+Q = 0, we can expect that 
in some limits the former can be approximated by one 
of the latter two. For the reflecting BC, with particles 
starting near +0 at time t — 0, the solution is 

p(M) = ^exp(-^-_), (A.6) 

where B is a constant. At a given time, the maximum 
value of this function is at x = and equals B; the 
x derivative at that point is zero. On the other hand, 
the maximum absolute value of the x derivative is at 
x ~ ^/DrI and is - ~ Byfrjgjt. The ratio 

between the maximum absolute value of the derivative 
and the maximum value of the function is U r/n/t. If this 
value is much larger than the ratio of the values of the 
derivative and the function at +0 given by Eq. (TOl) [this 
latter ratio is ~ (l/d)(ri R /r] L ) a }, then just a small shift in 
space of the solution of the reflecting problem, Eq. (IA.6I) , 
will be sufficient to satisfy the BC (|A.4|) . so in this situ- 
ation, which is observed when y/rjnjt ^ (l/d)(rjn/VL) a , 



or for 

t^t c = VR d 2 (^) ~^4~ 2a , (A.7) 

the solution of the reflecting problem, Eq. (|A.6[) ■ is an 
adequate approximation of the solution of the original 
problem with the BC (IA.4I) . On the other hand, consid- 
ering the solution of the absorbing problem, 

KM) = ^ex P (-^), (A.8) 

where C is a constant, the maximum value of the func- 
tion is ~ Cy/Dn/t ~ C /(t-y/rjii), the maximum value of 
the derivative is C/i 3 / 2 , and the ratio of the latter and 
the former is \J rfu/i, as in the reflecting case. But now, 
in order for the solution with the BC (|A.4[) to be ob- 
tainable by a small perturbation of that of the absorbing 
problem, we need this ratio to be much smaller than that 
specified by Eq. (|A.4|) . which happens for t t c , with t c 
still given by Eq. (|A.7[) . Thus, the solution of the Robin 
boundary- value problem (|A.5j) - (|A.4[) crosses over from 
being "reflccting-like" at t -C t c to being "absorbing- 
like" at t 3> t c . Note that the "reflecting- like" solu- 
tion corresponds to the particle number on R staying 
nearly constant, whereas the "absorbing-like" solution 
corresponds to this particle number decreasing rapidly 
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(as i -1 / 2 ) as the particles are "absorbed" by the inter- 
face, "transferred" to the left side, and "flow" to the wall. 

However, it should be remembered that our consider- 
ation is only valid in the interval t\ <C t <C t%. When 
a = 0, t c ~ t2, so the "reflecting-like" solution should be 
observed in the whole range of validity, although some 
"traces" of the "absorbing-like" solution can appear as 
one approaches £2- On the other hand, when a > 1/2, 
t c < ii, so the solution is "absorbing-like" in the whole 
interval. Of course, since for t < t\ the particles do not 
yet feel the walls, the distribution should be the same as 
in the reflecting case, that is, "reflecting-like" ; the tran- 
sition should occur around t\. It is only for < a < 1/2 
that the crossover should actually be observed in the in- 
terval ti < t < <2. If we choose t ~ y/hh, then for 
a = the particle distribution in the high-viscosity re- 
gion should resemble Eq. (IA.6I) ; for a — 0.5 and a = 1 it 
should be similar to Eq. (|A.8|) : finally, for a — 0.25 this 
time equals t c and the distribution should be a "hybrid" 
of these two cases. Figure [A~T1 showing the distributions 
for t]l = 2, T)r = 2000, t = 10 5 indeed matches these ex- 
pectations. Even though t — 10 5 is about 2.5 times larger 
than \Jt\t2 ~ 40000 (this larger time is chosen to have 
broader distributions, so their features are more visible), 
the crossover is broad enough that this time can still be 
considered to correspond to the crossover for a = 0.25. 
Note that the distributions in the low-viscosity region 
are linear, as given by Eq. (|A.1[) ; note also that a dif- 
ferent scale had to be used for the low-viscosity region, 
since the values of the particle density in that region are 
much lower than in the high-viscosity region. Thus, for 
any a, the fraction of particles in the low-viscosity re- 
gion is very small in this time range; contrast this with 
the reflecting case, where, for example, in the isothermal 
case there are more particles in the low-viscosity region 
than in the high- viscosity region at all times (except when 
t — > 00, these numbers become equal). To obtain the plot 
in Fig. IA.1| the variant of the MC I algorithm with the 
waiting time drawn from the distribution (1721) is used 
providing significant computational time savings for this 
large viscosity ratio. 

We are now in a position to obtain the FPT distribu- 
tion for the left wall, Ph{t)- In the "reflecting-like" sit- 
uation, considering the solution with the reflecting BC, 
Eq. (|A.6j) . by itself, the x derivative at +0 is, of course, 
zero. The next iteration is obtained by inserting the value 
of the function at this point, B/y/t, into the BC (|A.4|) . 
which gives ~ B(r] H /r] L ) a /(dVt). Using Eqs. (gBJl and 

4A3J, 



Px.it) 



BDr(vr/vlY 
ds/t 



drift 1 / 2 ' 



(A.9) 



Here we have used the fact that wl is always ~ 1. In the 
"absorbing-like" case, the x derivative at +0 is C/t 3 / 2 , 
and in the same way we obtain 




P L (t) 



C 



(A.10) 



FIG. A.l: (color online) The particle distribution functions 
for tjl — 2, r/R — 2000, and t — 10° with absorbing boundary 
conditions obtained using the variant of the MC I algorithm 
with the waiting time drawn from distribution (|72p using 10 6 
realizations. All particles start at the interface in the middle 
of the system. Different scales are used for the left and right 
parts of the plot: The left axis corresponds to the left part 
(with r/L = 2) and the right axis corresponds to the right part 
(with r/R = 2000). The distributions are normalized so that 
the area under the curves is unity, regardless of the fraction 
of particles left. 



Thus, the P(t) dependence should always be a power 
law; the exponent can be —1/2 in the "reflecting-like" 
situation or —3/2 in the "absorbing-like" situation. In 
Fig.[OJ we plot P L {t) for tj l = 2 and r] R = 2000. We see 
that indeed, for 0.5 < a < 1 (which corresponds to the 
"absorbing-like" situation at all times), the dependence 
is ~ i~ 3 / 2 . For a = (the "reflecting-like" case), the 
i -1 / 2 fit is adequate, although the slope does get steeper 
close to £2, as the crossover is starting to be "felt." For 
a smaller viscosity ratio, as in Fig. [71 the crossover is felt 
over the whole intermediate range of times, which is prob- 
ably what causes the apparent t~ 8 dependence in that 
case. Finally, for a = 1/4 the curve cannot be fitted by 
a simple power-law function; apparently, the crossover is 
so broad that the already large viscosity ratio needs to be 
increased by at least another order of magnitude before 
the distinct t~ x / 2 and t~ 3 / 2 regions can be observed. 

To obtain the prefactors (or the coefficients B and C) 
in Pi,(t), we assume that the number of particles on R, 
Nr, at the beginning of the intermediate time interval 
(t ~ ti) is still the same as at t — > 0, which, in its turn, 
is the same as in the reflecting case, since at t — > the 
particles have not reached the boundaries yet. Note that 
this does not mean that the ratio Nr/Nl at t ~ t\ is 
the same as at t — > 0, because particles on L do start 
to be absorbed around t\, so Nl does change by ~ t\. 
This means that while Eq. (IllOp can be used to get Nr 
at t ~ ti, Nl needs to be replaced with Nq — Nr, where 
Nq is the total initial number of particles on both sides. 
Two cases need to be considered. For a < 1/2, Eq. (|110[) 
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FIG. A.2: (color online) The FPT distributions for t}l = 2, 
r\B, = 2000, and d = 50 (dashed lines). All particles start at 
the interface in the middle of the system. The solid lines are 
the data for d = 100, rescaled appropriately (i.e., the time 
is divided by 4 and the value is multiplied by 4). The data 
are obtained using the variant of the MC I algorithm with 
the waiting time drawn from the distribution (|72p using 10 6 
realizations. 



(modified as above) gives Nr ~> N n — Nr, so 

N r /N ~ 1. (A.11) 

On the other hand, for a > 1/2, 

Nr/N ~ (rt R /VL) 1/2 - a « 1. (A.12) 

The first case is "reflecting-like" at t\ and until t c ; nor- 
malizing the particle distribution function so the area un- 
der the whole curve (L and R) is unity initially [and there- 
fore, given Eq. (|A.11|) . the area under R is still approx- 
imately unity at ~ t\], we obtain (B / 'v^i) V D Rt± ~ 1, 
thus B ~ 1/ \JDr ~ y/r]R, and using Eq. (|A.9|) , 



P L (t;a<l/2) 



a— 1/2 



a-1/2 



(A.13) 



After the crossover, Pz,(t) oc £ 3 / 2 ; by matching to 
Eq. (|A. 13|) at the crossover point t c , 



t a t 1/2 ~ c 

P L (t;a < 1/2) 



(A.14) 



The case a > 1/2 is "absorbing-like" for all intermedi- 
ate t; the area under the curve given by Eq. (|A.8[) at 
< = <i is ~ CDrJ^ 2 ~ C /(rjRt 1 ^ 2 ), so, according to 
Eq. (|A.12[) . the normalization condition is C/(?7ftiy 2 ) ~ 

(VR/ VL) 1/2 ~ a , thus C - f?B 2 ~ a *i /VL 2 ~ a > and > usin S 
Eq. (pOfj]) . 



Pi(t;a>l/2) 



fc l 6 2 



(A.15) 



Note that Eqs. (fA~l~3)) and (|A~1~5|) are not necessar- 
ily valid, even approximately, at t ks t\, mostly because 
around t±, and until time t* = (3ti, where /3 is a suffi- 
ciently large factor, the particle distribution on L may 
still be very different from that given by Eq. (IA.1[) . 
[Strictly speaking, this means that using t ~ t\ in consid- 
erations leading up to Eqs. (|A.13D — (|A.15|) is not entirely 
correct, but since for all practical purposes /3 can be cho- 
sen as a large, but finite and constant numerical factor, 
this does not affect the scalings, apart from at most log- 
arithmic corrections.] For a > 1/2, the number of parti- 
cles on R relative to the total initial number of particles 
in the system Nq is always <C 1. On the other hand, 
using Eq. (IA.3|) , the number of particles on L (again, rel- 
ative to A^o) at t = t* (at which the above considerations 
are already valid but still as close to ti as possible) is 



JVi(i*)~ A(f)d. 



P L {t*)d 2 



< 



"1/2 



< 1. 



(A.lf 



Thus by time t* most particles have already disappeared 
from the system. This absorption of particles must have 
happened around time ~ t\ over the time interval ~ t\. 
Thus, Pc(ti) - \jt\. However, Eq. (|A.15|) predicts 
Pl(*i) - (l/*i)(ii/*2) a_1/a « l/h . So, on top of the 
power law predicted by Eq. (|A.15I) there should be an 
additional peak (or "bulge") at ~ t\ of height ~ \jt\: 



P L (t 1 \a>l/2)~l/t 1 . 



(A.17) 



This "bulge", in fact, contains nearly all particles, al- 
though the rest of the distribution is still important when 
calculating the MFPT, as we will see. On the other 
hand, for a < 1/2 only a small fraction of all parti- 
cles, ~ (vr/vl) " 1 ' 2 , are on L initially. Even if all of 
them are absorbed over time ~ t\, the resulting FPT 



density is — {r]R/r]L) a ~ L/2 [h ~ *2 ~^ 2 Ai +1 ^ 2 , which is 
the same as what Eq. (|A.13|) at t = t\ gives, so there 
is no additional "bulge" . In fact, in Fig. IA.2I there ap- 
pear to be "bulges" for all a. However, it should be 
kept in mind that our discretization is not entirely ade- 
quate for this problem with a very high viscosity ratio: 
at t = ti = 1250, the width of the distribution on R 
is less than one, that is, less than one mesh step. In 
the same figure, we also show the results obtained for 
d = 100, with the resulting distributions appropriately 
rescaled. The "bulges" on the corresponding curves with 
a < 1/2 are smaller, so it is reasonable to assume that 
they largely disappear in the continuum limit. On the 
other hand, the "bulges" on curves with a > 1/2 are of 
the same size for d = 50 and d = 100 and therefore are 
"real" , rather than discretization artifacts. 

Using our results for FPT, we can now analyze the 
expression for the MFPT we have obtained [Eq. (|138[) ]. 
by decomposing the MFPT into parts associated with 
different regions of the FPT distribution. For a > 1/2, 
the weight of the bulge is nearly one, while the weight of 



which is the same expression as Eq. (|A.14|) . 
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the power-law region is 



P L (t)dt 



ia +a+l/2-a 
L 1 L 2 



dt 



a-l/2 



< 1. 



(A.18) 

On the other hand, the average time over the "bulge" is 
~ t\ , but the average time over the power-law region is 



which is much larger. The total MFPT is 



(A.19) 



(|) a_1/2 x(W 2 )~t 1 + t?4- c 



TL 



1 x h + — 



(A.20) 

where the first term comes from the "bulge" and the sec- 
ond one comes from the power-law region at intermediate 
times. Note that the second term is always larger than 
the first one (except for a = 1 , when they are equal) , and 
therefore it is the second term that determines the char- 
acter of the MFPT dependence, which indeed coincides 
with Eq. (|138l) . For a < 1/2, there are two power-law 
regions. The weight of the first, "reflecting-like" region 
is 



2 ,.. dt 
, tft 1 / 2 



r 



a— 1/2 ,1/2 



(A.21) 



The weight of the second, "absorbing-like" region is 



The contribution of the second term dominates (except 
for a = 0, when these contributions are equal, even 
though the second region "shrinks into a point"). Again, 
this second term coincides with Eq. (|138[) . Note that for 
both a < 1/2 and a > 1/2, contributions of the parts 
of the distribution with t > ti are significant. Particles 
spending such a long time before being absorbed neces- 
sarily spend most of that time in the high- viscosity re- 
gion. Note in this regard that the term in Eq. f|l 34[) asso- 
ciated with the time spent on R dominates for < a < 1 
when rjn/r/L — > oo, and even at a = 1 the terms cor- 
responding to the time spent on L and on R contribute 
equally (by order of magnitude). 

It is interesting to compare the particle distributions 
in Fig. IA.1I with those obtained for the same T)l, rjn, d, 
and t for reflecting boundaries. The crucial difference 
is that, while in the absorbing case the flow of parti- 
cles is from the high- viscosity to the low- viscosity region, 
the opposite is true in the reflecting case, which follows 
from the fact that for tjr/til > 1 and for all a, the ratio 
Nr/Nl grows with time [see Fig.[5fc)], yet Nl+Nr stays 
constant. As mentioned, this growth is proportional to 
yt; given that it occurs between t\ and and using the 
short-time [Eq. f|l 10[) ] and long-time [Eq. (jllip ] values of 
this ratio, we get 



Nr 
N l 



1/2- 



t 



-t 



1/2 



(A.26) 



For a < 1/2 this ratio is much larger than unity at all 
times, and the fraction of particles on R is 



+a j-1/2 
l l l 2 



£3/2 



-dt 



.iaj.l/2-a 



1/2 



(A.22) 



So both regions carry equal weight. In fact, the largest 
contribution comes from the region around t c — the same 
result is obtained if the lower limit of the first integral 
and the upper limit of the second integral are replaced by 
arbitrary values, as long as the first one is much smaller 
than t c and the second one is much larger than t c . To 
put it differently, of all intervals [t a , f ft] with equal ti,/t a , 
the one containing t c will contribute the most. How- 
ever, other regions cannot be ignored when calculating 
the MFPT. The average time over the "reflecting-like" 
region is 



f, to t^dt 



m ,2a,.l-2a 
i c ~ lj l 2 , 



r t- i / 2 dt 

while that over the "absorbing-like" region is 
J t t2 t- x l 2 dt 



(A.23) 



f t t2 t~ 3 / 2 dt 



+1/2A/2 +a+l-a 
L c L 2 L l l 2 



The resulting MFPT is 

t l ~ 1 X t\ a t\- 2a + 1 x t^t\- 



(A.24) 



(A.25) 



N L +N, 



1 - 



R 



Nl 
Nr 



1 - 



t 



1/2- 



-t 



-1/2 



(A.27) 



The derivative of this is equal to the flux of particles into 
the right region, which in its turn is proportional to the 
space derivative of the particle distribution function at 
x = +0. Thus, 



1 dp(x,t) 



Ox 



t 



x=+0 



1/2-a 



dp(x,t) 



dx 



t 



x=+0 



t. 



1/2-a 



VRt 



-3/2 



-3/2 



(A.28) 



(A.29) 



Assuming that the distribution on R can be approxi- 
mated by one half of a Gaussian curve with area ~ 1 
[Eq. (|A.6p with B ~ \ j \JDr ~ y 7 ?/]?], the maximum ab- 
solute value of the x derivative is ~ ?7_r/£- This is much 
larger than the absolute value of the derivative given by 
Eq. (|A.29[) for any t ^> ti. Thus, arguing as in the absorb- 
ing case, the Gaussian with the peak at zero [Eq. (|A.6[) ] 
is indeed a good approximation in this case. 

Consider now a > 1/2. In this case, two situations are 
possible: for 



t > t c 



-2a+2a- 



(A.30) 
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Nr/Nl is still large, and the considerations above for 
a < 1/2 mostly apply. In particular, Eq. ([A. 291) and 
the expression for the maximum absolute value of the x 
derivative of the trial Gaussian are still valid, and the 
absolute value of the former is much smaller than the 
latter for any t 3> t cr , so the Gaussian with a peak at zero 
still approximates the distribution well. On the other 
hand, for t < t cr , N R /N L < 1. This means that N L 
stays nearly constant and thus p(x, t) on L, besides being 
nearly constant in space, is also nearly constant in time 
and 



/ x 2 1 
p(x, t) w - ~ ~, x < 0, t <C t cr . 
a a 

Then, according to the jump condition (j4"4")l , 



p(x 



-0,t)~* 

d V Vl 



(A.31) 
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and is likewise time- independent. The solution of the dif- 
fusion equation satisfying the condition p(+0, t) = const 
for t > is the complementary error function erfc [64j . 
and we get 



p(x, t) ~ — a /—erfc 



x>0. 



(A.33) 



Note that Nr ~ vt, as expected, and also that the in- 
flection point is now at x = +0 (indeed, dp/dt = at 
a; = +0 implies d 2 p/dx 2 = 0). 

Similar to the absorbing case, we can consider the par- 
ticle distributions at t ~ y/t\t2- We expect the Gaussian 
distribution for a between and 0.5, the complementary 
error function with the inflection point at for a = 1 
and a "hybrid" situation for a = 0.75. In Fig. IA.3| we 
show the results obtained using MC I simulations with 
t\ w given by Eq. (|73p . These results confirm our expecta- 
tions. The a = 0.75 curve resembles the a = 1 curve, but 
in the former case the inflection point is clearly shifted 
to the right of the interface. 

Comparing the results in Figs. IA.1I and IA.31 it is in- 
teresting to note a kind of "mirror symmetry" between 
them. For example, in the absorbing case the distribu- 
tions with a between 0.5 and 1 are very similar in the 
high- viscosity region, and in the reflecting case the distri- 
butions with a between and 0.5 have this property. In 
the absorbing case, a — 0.25 is "hybrid," while in the re- 
flecting case, a = 0.75 is "hybrid." Finally, in the absorb- 
ing case the distribution for a = 1 is given by Eq. (IA.8I) , 
which is the derivative of the Gaussian distribution, and 
in the reflecting case the distribution for a = 1 is given 
by the complementary error function [Eq. (|A.33|) ]. which 
is the integral of the Gaussian distribution. 



Appendix B: Proportion of the number of steps on 
left and right for a random walk absorbed at a 
particular wall 

Using the MC I algorithm with all move attempts suc- 
cessful (Sec. IIV A 2[) reduces the problem of particle dif- 



Q = 

a = 0.50 
a = 0.75 

o = l 




FIG. A. 3: (color online) The particle distribution functions 
for t]l = 2, rjR = 2000, t — 10 5 with reflecting boundary 
conditions obtained using the variant of the MC I algorithm 
with the waiting time given by Eq. (|73[) using 10 6 realizations. 
All particles start at the interface in the middle of the system. 



fusion in a system with a diffusivity interface to a lattice 
random walk that is unbiased everywhere, except per- 
haps at the interface site. Depending on what side of 
the interface the particle is on, a single step of this ran- 
dom walk takes a different amount of time, but this can 
be taken care of separately, by working in terms of the 
numbers of steps on each side and then multiplying by 
the respective time steps. In this way, our problem is 
reduced to that of an ordinary lattice random walk with 
a constant diffusion rate that starts at the site in the 
middle between two absorbing walls and can be biased 
at that site, but nowhere else. We first consider the case 
with no bias even at the middle site (corresponding to 
the Ito calculus) and then generalize. 

First, consider the diffusion on a lattice with mesh step 
a = 1 of an ensemble of particles released half-way be- 
tween walls at x — ±d/2, with d/2 integer and the initial 
position of each particle given by xq = 0. By symmetry, 
the average mean first passage times (MFPT) to either 
wall (the conditional MFPTs) are equal and thus equal 
to the MFPT to either wall (the unconditional MFPT). 
If the time is measured in terms of the number of steps, 
then the MFPT is the average number of steps, which is 

M 



<p_ 

4 ' 



(B.l) 



To proceed, we can decompose the full trajectory of 
a particle into segments separated by its visits to the 
middle site. Each such segment is entirely either in the 
left half or in the right half and for brevity, we refer 
to these groups of segments as left and right segments. 
Consider the ensemble of trajectories ending on the right 
wall. Taking one such trajectory and replacing all of the 
right segments (except for the last one, of course) with 
its mirror-image left segment and vice- versa produces an- 
other valid trajectory belonging to the ensemble. There- 
fore, up to the last segment, there should be an equal 
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number of left and right segments on average and they 
should have the same average duration. Thus, we can di- 
vide the trajectory into two parts. The first part contains 
all but the last segment and on average, the time spent 
in both halves of the interval is the same, but the second 
part, consisting of only the last segment, is entirely in 
the right half. 

The average time corresponding to the second part of 
the trajectory can be obtained by recasting this problem 
as an equivalent one. Consider the interval consisting of 
just the right half of the original interval. If both bound- 
aries of this interval are treated as absorbing, the particle 
is released right next to the boundary corresponding to 
the middle site and the conditional MFPT, given that the 
particle is absorbed by the other boundary, is calculated, 
then all trajectories which return to the middle site will 
be discarded. While the probability of being absorbed by 
a wall vanishes in the limit when the distance from the 
point where the particle is released to the opposite wall 
approaches zero (on a lattice this probability is ~ a/d), 
the MFPT to the far wall is well-defined and is, for the 
continuum problem [65j | . 



and the average number of steps on the opposite side is 



''opp 



2 U1 = IT 



(B.6) 



Consider now the case when there is bias at the inter- 
face, so a particle at the middle site chooses to move left 
and thus start a left segment with probability or to 
move right and thus start a right segment with probabil- 
ity Pi+ = 1 -Pi_- These probabilities for MC I are given 
by Eqs. (f?4"]) and ([75]) . Each of these segments is unbiased 
and thus, statistically speaking, unmodified compared to 
the case of no bias. In particular, the distribution of the 
number of steps in each segment is the same and the 
probability to reach the wall is the same as well. This 
means that Eqs. (jB.ll) . (|B.3I) . and (|B.4|) are still valid. 
The only thing that changes is the fractions of left and 
right segments in the first part of the trajectory, which 
are now p'®_ and p'i + , respectively. Then, for trajectories 
ending on the right wall, the average number of steps on 
R is 



T-2 
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where I = d/2 is the distance between the walls (the 
length of the half-interval under consideration). Using 
Eq. ([55]) with p s = 0, a — 1 and r = 1, we convert this 
result into that for the average number of steps in the 
corresponding lattice problem, 



112 = 12- 



(B.3) 



This is the average duration of the second part of the 
trajectory. The average duration of the first part is 



flRR 



and on L, 
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Similarly, for trajectories ending on the left wall, the av- 
erage number of steps on L is 



ni = n 



no = — . 
6 



(B.4) 



n L L = pf-ni + n 2 = — (1 + 2p' x °_), 



(B.9) 



Given that on average, the first part contains an equal 
number of steps on L and on R, but the second part 
is entirely on the same side as the wall of interest, the 
average number of steps on the same side as the wall is 



1 

2 ni 



712 = ¥' 



(B.5) 



while on R, 



n L R = Pi+ni = — p'i + 



(B.10) 
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